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INTRODUCTION 


Heat  exchange  between  the  human  and  the  environment  has  always  been  a  topic  of 
great  interest  as  this  is  one  of  the  essential  manifestations  of  homeothermy.  This  interest  has 
intensified  over  the  past  few  decades  as  the  understanding  of  the  mechanisms  involved  are 
spurred  on  by  human's  venturing  into  extreme  environments,  e.g.,  outer  space.  Even  less 
extreme  environments  may  pose  life  threatening  challenges  to  humans  and  the  literature  is 
laden  with  examples  of  both  heat,  [1]  and  cold,  [2]  related  casualties. 

Physiological  studies,  in  which  human  subjects  are  exposed  to  extreme  environmental 
conditions  are  essential  for  collecting  data  on  the  actual  thermal  behavior  of  men  and 
women.  These  studies  are  employed  to  generate  detailed  and  reliable  databases  and  to  test 
thermoregulation  theories.  Much  information  may  also  be  gained  by  formulating  models  which 
simulate  qualitatively  the  thermal  behavior  of  the  human  body.  The  chief  advantage  of  these 
models  is  their  ability  to  predict  and  point  out  trends  and  limitations  while  avoiding  the 
dangers  and  cost  involved  in  actually  performing,  time  consuming  and  sometimes  hazardous 
experiments.  Their  inherent  disadvantage  resides  in  the  necessity  to  verify  their  predictions. 
A  variety  of  models  simulating  human  thermal  behavior  have  been  developed.  These  range 
from  models  of  single  organs,  e.g.,  [3-5]  to  models  of  the  entire  body,  e.g.,  [6-9]. 

In  this  report  a  detailed  model  of  an  extremity  exposed  to  cold  weather  is  developed. 
The  reasons  for  choosing  an  extremity  are  two-fold:  (a)  a  model  of  an  extremity  may  serve 
as  a  "building  block"  for  other  elements,  and,  (b)  the  extremities  are  usually  the  most 
vulnerable  body  elements  particularly  in  cases  involving  operations  in  cold  weather. 

The  extremity  is  depicted  as  a  right-angle  cylinder  in  which  heat  flows  in  both  the 
radial  and  axial  directions.  Around  the  entire  external  surface  of  the  cylinder  different  layers 
of  insulation  may  be  applied  through  which  heat  is  exchanged  with  the  environment.  Heat  is 
also  exchanged  internally  by  conduction  and  with  the  blood  flowing  both  in  the  major  blood 
vessels  and  in  the  vessels  of  the  capillary  bed  Counter-current  heat  exchange  between  the 
major  blood  vessels  is  also  taken  into  account. 

The  model  is  presented  as  a  consistent  set  of  energy  balance  equations  and  is  solved 
by  a  finite-difference,  alternating  directions  numerical  scheme  employing  the  Thomas 
algorithm.  This  scheme  has  been  tested  extensively  for  stability,  convergence  and  accuracy. 
It  was  also  run  for  a  number  of  cases  to  demonstrate  its  fundamental  capabilities. 
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ANALYSIS 

Energy  balance  in  a  right  -  angle  circular  cylinder  depicted  in  Fig.  1  is  expressed  by: 


pc 


dT' 

at * 


=  k 


_J_  (r.  ar;,  +  I 

3r  *  3r'  a 


v*  Cb  ( r;  -  r* )  +  U.  ( I*;  -  r* )  +  u„  ( r;  -  r* ) 


(1) 


where  the  term  on  the  left  hand-side  represents  the  rate  of  change  of  stored  energy  and  the 
terms  on  the  right  hand-side  express  radial  and  axial  heat  conduction,  metabolic  heat 
generation  rate,  heat  exchange  with  the  capillary  bed  and  heat  exchange  with  a  large  artery 
and  a  large  vein,  respectively.  All  properties  and  variables  are  defined  in  the  Glossary  and 
asterisks  indicate  dimensional  variables. 

The  following  boundary  and  initial  conditions  are  specified  for  the  problem: 

At  the  center  of  the  cylinder  an  adiabatic  condition  is  formulated  to  satisfy  symmetry 
requirements: 


—  =  0  9  r*  =  0  (2) 

dr * 


On  the  outer  circumferential  surface  of  the  cylinder  heat  is  exchanged  with  the  environment 
by  convection: 


(  *0*  -  T*  ) 

or  * 


9  r*  =  R 


At  the  base  of  the  cylinder  a  variable  temperature  is  assumed: 


(3) 
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*1  ( f) 


9  z*  =  0 


(4) 
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At  the  tip  of  the  cylinder  convective  heat  exchange  with  the  environment  occurs: 

=  ~  (r0*  -  t*)  •  «*  =  l  (5) 

a**  k 

Initially  the  temperature  distribution  in  the  cylinder  is  expressed  by  an  arbitrary  function: 
IP*  =  Tl  (r* ,  z*)  9  t*  =  0  (6) 


Equation  (1)  contains  terms  expressing  the  two  modes  by  which  heat  is  exchanged 
between  the  tissue  and  the  circulatory  system.  In  these  expressions  T’  and  T*  represent 
arterial  and  venous  temperature  distributions,  respectively.  It  is  assumed  that  the  cylinder  is 
traversed  by  one  each  of  these  large  blood  vessels  in  the  axial  direction  (Fig.  1).  In  addition 
to  exchanging  heat  with  the  surrounding  tissue,  these  two  vessels  are  also  coupled  by 
counter-current  heat  exchange.  Two  separate,  but  coupled  heat  balance  equations  are  now 
written  for  the  large  artery  and  the  large  vein,  respectively: 

inCb^m,  in  “  out  Cb^*,  out  + 

(7) 

Ju.  (T*  -  r; )  dv  +  hmv  ( iv*  -  r;>  -  / wbcbTl dv 


dT; 

~dP 


and, 


M„Ct 


dT* 

df 


K,  in  Cb^v,  in  1K’,  out  Cb  out  + 


(8) 


J"uv(T*  -  T*)  dv  +  hay,  (T*  -  T*)  +  fwbcbT*  dv 


The  terms  on  the  left  hand-side  of  Equations  (7)  and  (8)  represent  the  rate  of  change 
of  the  average  amount  of  energy  stored  in  the  blood  contained  at  any  instance  in  these  two 
vessels.  The  terms  on  the  right  hand-side  of  these  equations  represent,  respectively,  the 
enthalpy  flux  into  and  out  of  the  control  volume,  the  contribution  to  the  heat  balance  due  to 
the  heat  exchange  with  the  average  temperature  of  the  surrounding  tissue,  T,  and  the 
contribution  due  to  counter-current  heat  exchange  between  the  two  large  blood  vessels.  The 
remaining  terms  in  Equations  (7)  and  (8)  indicate,  respectively,  the  drainage  into  and  the 
collection  from  the  tissue  by  capillary  perfusion  as  the  two  vessels  traverse  the  cylinder.  The 
heat  transfer  coefficients  between  the  blood  vessels  and  the  surrounding  tissue  (ua  and  uv), 
and  between  the  two  blood  vessels  (h^)  are  derived  in  Appendix  B. 
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Two  additional  mass  conservation  equations  are  required  for  both  large  blood  vessels: 


K.in  =  K.out  *  f  <9) 

and, 


K 


in 


(10) 


Prior  to  applying  a  numerical  solution  to  the  coupled  Equations  (1 ),  (7)  and  (8),  subject 
to  boundary  and  initial  conditions  (2)-(6),  these  equations  are  rewritten  in  dimensionless 
forms: 


dr  =  JL  A  JL  (r  — )  +  —  d*T 

dx  ab  r  dr  '  dr  1  a2  dz 2 
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/"■ 
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R3 
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[  K,  Id  ^b  ®V,  in 
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ar 

dr 


=  0 
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uv(T  -  2V)  dv  +  hav 


@  r  =  0 


(11) 
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T  =  Tx  (  t )  9  z  =  0 

^  =  Bix  (r0  -  T)  9  z  =  1 

dz 

r=ri(r/*)  «t  =  o 


where, 
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r J 
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(16) 
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(26) 
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(27) 
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Jfc 


(28) 


and, 

Bi,  =  ^2-^  (29) 

1  it 

A  finite-difference  solution  is  formulated  for  the  above  set  of  dimensionless  equations. 
The  cylinder  is  divided  into  four  radial  compartments  depicting  the  core,  muscle,  fat  and  skin, 
respectively,  as  shown  in  Fig.  2.  Each  of  these  compartments,  the  boundaries  of  which  are 
determined  by  anatomical  and  anthropometric  considerations,  may  be  further  subdivided 
radially  according  to  the  required  details  of  the  temperature  variations  in  the  cylinder.  Axial 
divisions  are  uniformly  spaced.  A  cross  section  of  a  typical  control  volume  is  shown  in  Fig. 
3. 

As  a  first  step  in  the  numerical  solution  ,  Equation  (11)  is  multiplied  by  a  hollow 
cylindrical  volume  element  of  thickness  dr  and  length  dz  and  integrated: 


/  £*drd* 


_a_  1  Jl  ,  dT)  +  _l_  *T  + 
ab  r  dr  dr 1  aa  dz 2 


(30) 


▼  •  [  q  +  ( rr  +  um)  ( ra  -  t  )  +  ( rT  -  t  )  ]  j  rdrdz 

The  temporal  derivative  on  the  left  hand-side  of  Equation  (30)  is  calculated  by  a 
forward  difference.  In  evaluating  this  derivative,  a  half  time  step  is  assumed  to  facilitate  a 
solution  of  this  two-dimensional  problem  by  the  method  of  alternating  directions  [10].  In  this 
method  the  solution  of  the  resulting  set  of  difference  equations  is  performed  in  two  half  time 
steps:  first  the  temperatures  in  one  direction,  e.g.,  radial,  are  calculated  for  the  first  half  time 
step,  based  on  the  values  of  the  temperatures  at  the  current  time  in  the  other  direction.  Next, 
the  temperatures  at  the  other  direction,  e.g.,  axial,  are  evaluated  for  the  next  half  time  step 
based  on  the  values  obtained  for  the  first  spatial  direction  in  the  first  half  time  step.  This 
completes  a  full  time  step,  and  these  values  are  used  to  initiate  the  next  full  time-step 
iteration. 
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Figure  3:  Cross  section  of  a  typical  control  volume  of  the  cylindrical  model 
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Accordingly,  the  integration  of  Equation  (30)  yields: 
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where  superscript  n  indicates  full  time  steps  and  i  and  j  indicate  radial  and  axial  divisions, 
respectively. 

Canceling  identical  factors,  rearranging  and  redefining  the  temporal  and  spatial 


divisions  by: 

M 

< 

ll 

«r 

(32) 

hM  =  A  z 

(33) 

»  ■  % 

(34) 

yields: 
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Equation  (35)  is  the  general  finite-difference,  or  discretized,  equation  for  the  tissue 
temperatures.  This  equation  may  more  conveniently  be  rewritten  in  the  following  form: 


1  +  ft 


_a_  _2_  +  ▼  (»r  ♦  P«  +  uT) 

ab  hi  2 


For  simplicity  a  certain  notation  convention  is  adopted  in  Equation  (36)  regarding  the 
spatial  indices.  Accordingly,  whenever  a  nominal  spatial  index  i  or  j  occurs,  it  is  dropped  out 
from  the  equation.  This  leaves  only  stepped  indices  to  be  specifically  indicated,  e.g.,  Ti  j+1  is 
represented  by  V  etc. 

Equation  (36)  is  evaluated  for  all  nodal  points  in  the  cylindrical  domain.  The  process 
involves  substitution  of  the  boundary  conditions  (equations  (14)-(i7)},  accounting  for  the 
dissimilar  nodal  spacings  at  the  boundaries  between  the  various  tissue  compartments  in  the 
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radial  direction,  substitution  of  zero  values  for  perfusion  near  the  external  boundaries,  etc. 
Details  of  the  derivation  are  presented  in  Appendix  A. 


This  process  yields  a  set  of  algebraic  equations  for  all  the  tissue  nodal  points.  Each 
equation  in  this  set  usually  includes  three  terms  for  the  radial  direction  and  three  additional 
terms  for  the  axial  direction.  An  additional  term  not  containing  unknown  tissue  temperatures, 
is  also  included  in  each  equation.  At  the  boundaries  in  both  directions  only  two  terms  are 
present,  yielding  tri-diagonal  matrices  for  the  algebraic  set  of  equations.  This  property  of  the 
set  of  equations  renders  it  solvable  by  the  Thomas  algorithm  [11].  To  apply  this  algorithm 
Equation  (36)  is  now  written  in  matrix  notation: 

[  J  -  6  ‘Ar  ]  •  |  Ti*  2  |  =  [  I  +  6  -As  ]  •  {  T*  }  +  6  *  |  3d*1  J  (37> 


where  [I]  is  a  unit  matrix,  A,  and  \  are  the  elements  of  the  tri-diagonal  matrices  in  the 
radial  and  axial  direction,  respectively,  and  {Sn+1/2}  is  a  one  dimensional  vector  containing  all 
remaining  terms  in  Equation  (35)  which  do  not  multiply  the  tissue  temperatures.  Derivation 
of  these  quantities  for  all  nodal  points  is  presented  in  Appendix  A  and  a  summary  of  all 
coefficients  is  given  in  Tables  1-3. 

It  is  noted  that  Equation  (36)  indicates  the  calculation  of  the  first  half  time-step  only. 
To  complete  a  full  time-step,  an  equation  similar  to  Equation  (36)  is  required: 

[  I  -  5  \A,  ]•  |  T/*1  J  =  [  j  +  6  -Ax  ]  .  |  T?*  |  +  6  •  |  <38> 


The  particular  formulation  employed  in  the  present  analysis  assumes  that  the  S-vector 
in  Equations  (37)  and  (38)  is  evaluated  only  once  per  full  time-step,  i.e.,  at  the  one  half  time- 
step.  It  is  then  maintained  constant  for  the  two  calculation  passes  in  both  radial  and  axial 
directions.  The  S-vector  contains  the  temperatures  of  the  large  blood  vessels  at  the  axial 
nodal  points  which  provide  the  thermal  coupling  between  the  tissue  and  the  circulatory 
system.  To  calculate  these  temperatures,  a  forward-difference  approximation  for  the  time 
derivatives  in  Equations  (12)  and  (13)  is  employed.  Two  additional  assumptions  are  made, 
regarding  the  average  arterial  and  venous  blood  temperatures  and  average  tissue 
temperature  appearing  in  these  equations.  First  the  average  temperature  of  the  blood  in  any 
one  of  the  two  vessels  is  taken  as  an  arithmetic  mean  of  the  two  enclosing  axial  nodal  points: 

?»-■§  ('fcj.  ♦  **,«*>  (39) 
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Table  2:  Coefficients  in  the  axial-direction  [AJ  matrix. 
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Similarly,  the  average  tissue  temperature  for  exchanging  heat  with  any  of  the  large  blood 
vessels,  is  taken  as  the  arithmetic  mean  of  the  tissue  temperatures  at  the  two  enclosing  axial 
nodal  points: 


T  =  |[  +  T(i,j-i)  ] 


(40) 


With  these  assumptions  the  integrals  indicated  in  Equations  (12)  and  (13)  are  replaced  by 
numerical  summations  to  yield: 


Ta  2  (j)  =  Ax|  [  -  2  Bm  ( j-x )  +  Am  mSOMl  -  YM  -  Hmv  T°  ( j )  + 

[  2  Bm  (J- 1  )  -  Aa  'SDMI  ]  Ta  (j- 1  )  +  Ya  -SXJM2  +  Hav -T*  ( j  )  J 


(41) 


and, 


rv  2  (j)  =  At  |  ~  +  2  -Br  (  j  -1  )  -2  ’Ay  ■  SDMI  -  Yv  -  H„ 


[  -2  -B„  (  j  -1  )  +  A^  ‘SDMI  ]  T?  (  j  -  1  )  +  -|  Yy’80M2  +  H^-T?  (j)  +  (42> 

m* 

±Ay-SOM3  | 


where: 
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(43) 
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(46) 


jr-i 


SCBfl  =  g  *U.j-l)’{rln  -  r?) 


(47) 


SDM2  = 


g  {  T(i,j-1)  *  T[l,j)  }-(r/+i  -r/) 


(48) 


and, 


ir-i 


SUM3  = 


g  W  ( i  #  )  *{  T  (  i ,  j-1 )  +  T(i,j)  }  •  (  r/+1  -  r/  ) 


(49) 


According  to  the  present  formulation,  the  symmetry  condition  at  the  centerline  of  the 
cylinder,  Equation  (2),  is  satisfied  by  excluding  the  first  one-half  division  in  the  radial 
direction,  Fig.  A.2.  However,  in  performing  the  summations  indicated  in  Equations  (47)-(49), 
the  contribution  of  this  region  is  included  in  order  that  mass  conservation  requirements  be 
satisfied. 


PHYSICAL  AND  PHYSIOLOGICAL  PARAMETERS 


Three  interrelated  groups  of  parameters  are  required  for  calculating  the  temperature 
field  in  the  model.  These  are:  (a)  geometrical  parameters,  depicting  the  anatomical  details 
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of  the  modeled  organ,  (b)  thermophysical  parameters,  representing,  primarily,  the  transport 
properties  of  the  tissue,  and,  (c)  physiological  parameters  simulating  variables  such  as  blood 
flow  or  metabolic  heat  generation  rate. 

Accurate  and  detailed  information  on  these  parameters  is  not  available.  Moreover, 
individual  variabilities  among  humans  make  K  almost  impossible  to  formulate  a  universal  set 
of  parameters  for  the  model.  Nevertheless,  a  reasonably  accurate  set  of  parameters  may  be 
identified  for  the  purpose  of  studying  the  behavior  of  the  model. 

Table  4  lists  the  properties  used  in  the  model.  Data  are  given  for  the  four 
compartments,  or  organs  which  make  up  the  model,  i.e.,  core,  muscle,  fat  and  skin. 
Additional  data  are  given  for  blood.  Most  of  the  entries  in  Table  4  were  compiled  from 
References  [8]  and  [1 2].  Blood  perfusion  and  metabolic  heat  generation  rates  were  estimated 
as  follows.  According  to  Burton  [13],  the  average  blood  flow  in  the  finger  of  a  subject  "who 
is  comfortable  as  regards  the  temperature  of  the  surroundings"  is  in  the  range  of  15-40 
cc/min/100  cc  tissue.  We  assumed  the  lower  limit  of  this  range  to  be  representative  of  the 
basal  blood  flow  rate  in  the  unperturbed  finger.  Converted  into  SI  units,  and  assuming  Table 
4  value  for  blood  density,  this  basal  value  is  given  by  2.65  kg/m3  sec.  This  basal  rate  was 
used  for  calculating  the  organ-specific  values  by  assuming  the  geometrical  values  of  Table 
4  and  accounting  for  the  absence  of  blood  flow  in  the  fat  layer. 

Also  listed  in  Table  4  are  the  values  for  "nutritional"  blood  flow  rates  in  the  various 
organs  of  the  finger.  These  values  represent  the  flow  rates  in  a  fully  constricted  finger 
exposed  to  a  cold  environment.  Values  in  the  literature  for  this  condition  are  in  the  range  of 
0.3-1 .0  cc/min/100  cc  tissue  [13-15].  For  most  of  this  study  we  assumed  a  nutritional  blood 
flow  value  of  0.5  cc/min/100  cc  tissue  or,  0.0883  kg/m3  sec.  Values  for  the  different  organs 
cf  the  finger  are  listed  in  Table  4. 

Yet  another  set  of  values  relates  to  the  basal  metabolic  heat  generation  rate  in  the 
various  organs.  These  were  estimated  by  assuming  that  the  nutritional  blood  flow  rate  is 
maintained  for  the  purpose  of  supporting  the  metabolic  activities  of  the  tissues  under  all 
conditions.  Thus,  average  oxygen  extraction  rates  may  be  assumed  for  estimating  the  basal 
metabolic  heat  generation  rates.  According  to  Cooney  [16],  typical  oxygen  concentration 
levels  in  the  blood  are  0.195  and  0.145  liter  Oj/liter  blood  for  tissue  inlet  and  outlet, 
respectively.  The  average  caloric  equivalent  of  1  liter  of  oxygen  may  be  taken  at  20.9  kJ  (5 
kcal)  to  yield  the  basal  metabolic  heat  generation  values  listed  in  Table  4. 

In  listing  these  values,  one  adjustment  was  made  in  regard  to  the  metabolic  rate  of 
the  fat  layer.  Although  basal  blood  flow  rate  to  this  organ  was  assumed  to  be  practically  zero, 
some  metabolic  activity  could  still  be  assumed  for  this  organ.  Accordingly,  a  very  low  level 
of  5  W/m3  was  arbitrarily  assigned  to  the  fat  layer. 
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THERMAL 

CONDUC¬ 

TIVITY, 

SPECIFIC 

HEAT, 

DENSITY. 

W/m°C 

kJ/kg°C 

kg/m3 

1.064 

2.102 

1401 

0.418 

3.136 

1057 

0.204 

2.520 

900 

0.293 

3.780 

1057 

0.450 

3.899 

1060 

METABOLIC 

RATE, 


TO 

FINGER 

RADIUS, 


BASAL 

BLOOD 

FLOW 

RATE, 


NUTRI¬ 

TIONAL 

BLOOD 

FLOW 

RATE, 


kg/m3sec  kg/m3sec 


CORE 


MUSCLE 


170.5  5.195  0.173 


631.9  19.225  0.641 


247.4  7.526  0.251 


BLOOD 


Table  4:  Property  values  used  in  the  numerical  computations. 
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Heat  transfer  coefficients  used  to  represent  the  conditions  at  the  surface  of  the  finger 
are  listed  in  Table  5.  Four  combinations  are  considered:  bare  and  gloved  finger  in  either  still 
air  (free  convection)  or  windy  air  (15  km/h r).  The  values  were  calculated  by  standard 
engineering  equations  [17]  for  a  0.08  m  long,  0.015  m  diameter  cylinder.  A  distinction  was 
made  between  the  cylindrical  surface  of  the  finger  along  its  axis  versus  the  spherical-like  tip. 
The  glove  was  represented  by  a  3-layer  ensemble  depicting  a  2.86  mm  (0.09  in)  wool  layer, 
1  mm  of  still  air  gap  and  a  1 .27  mm  (0.05  in)  leather  shell.  Also  shown  in  this  table  are  the 
equivalent  clo  values  of  the  various  entries  which  conform  well  to  the  range  of  values 
measured  on  a  variety  of  gloves  [18,19]. 


TESTING  OF  THE  NUMERICAL  CODE 


A  rigorous  series  of  benchmark  tests  was  devised  and  carried  out  to  verify  the  stability 
and  convergence  of  the  numerical  code  written  for  the  model.  Programming  was  done  in 
Turbo  Pascal  Version  6.0  for  IBM-compatible  personal  computers.  Appendix  C  contains  a 
complete  listing  of  the  source  code  and  the  operating  instructions  for  the  program. 

In  the  first  group  of  tests,  all  physiological  parameters,  i.e.,  qm,  wb,  ua,  uv  and  hav  were 
set  to  zero.  This  rendered  the  problem  a  simple,  two  dimensional  heat  transfer  problem.  In 
tests  #1  -3,  the  heat  transfer  coefficients  on  the  surfaces  of  the  cylinder,  h,  and  hc  were  also 
set  to  zero  thereby  creating  an  adiabatic  cylinder  except  for  the  base  (z=0).  In  test  #1  initial 
and  boundary  temperatures  at  z=0  were  set  to  30°  C  and  the  program  was  run  for  200,000 
time  steps,  0.1  second  each.  Throughout  the  test  no  temperature  changes  were  observed 
anywhere  in  the  mesh,  as  is  to  be  expected.  In  tests  #2  and  3,  a  change  was  made  in  the 
boundary  condition  at  z=0  after  the  initial  100  time  steps.  In  test  #2,  run  for  100,00  time 
steps,  1  second  each,  the  change  was  from  20°  C  to  30°  C.  The  inverse  change  was  made 
in  test  #3  which  was  run  for  a  total  of  300,000  time  steps.  In  both  cases  the  temperatures 
anywhere  in  the  mesh  approached  the  boundary  temperatures  and  remained  stable. 

In  test  #4  an  active  heating  source  was  introduced  into  the  cylinder,  i.e.,  qm>0,  while 
still  maintaining  the  other  parameters  inactive  as  above.  Values  used  for  the  heating  source 
were  those  listed  in  Table  4  for  the  basal  metabolic  heat  generation  rate.  The  program  was 
run  for  a  total  of  2,000,000  time  steps,  0.1  second  each.  Mesh  temperatures  have  stabilized 
after  600,000  time  steps. 

The  results  of  the  cylinder  model,  with  an  internal  heating  source,  were  compared  to 
those  calculated  by  a  one-dimensional  analytical  solution  [20].  The  comparison  between  the 
surface  (external)  temperatures  as  calculated  by  this  model  and  those  of  the  analytical  model 
is  shown  in  Fig.  4.  It  is  noted  that  precise  comparison  of  these  two  cases 
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hc 

W 

h , 

W 

mz  °C 

m2  °C 

BARE  FINGER, 

STILL  AIR 

7.02 

(0.913  clo) 

9.46 

(0.672  clo) 

BARE  FINGER, 

15  km/hr  WIND 

67.2 

(0.095  clo) 

90.3 

(0.071  clo) 

GLOVE, 

STILL  AIR 

5.17 

(1.240  clo) 

8.09 

(0.792  clo) 

GLOVE, 

15  km/hr  WIND 

10.02 
(0.640  clo) 

13.09 
(0.490  clo) 

Table  5:  Heat  transfer  coefficients  for  combinations  of  wind  conditions  for  a  bare  and 
gloved  finger. 
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ANALYTICAL 

NUMERICAL 


DIMENSIONLESS  DISTANCE  ALONG  CYLINDER 


Comparison  of  temperature  distributions  along  the  finger  model  calculated 
by  an  analytical  model  [20]  and  the  present  model. 


is  not  possible  since  no  radial  temperature  variations  are  included  in  the  analytical  model. 
Nevertheless,  the  two  sets  of  plotted  results  seem  to  agree  well,  with  the  numerical  results 
slightly  under  predicting  the  analytical  ones  at  the  shorter  times  into  the  run. 

In  test  #5,  qm  was  reset  to  zero  and  h,  and  hc,  the  heat  transfer  coefficients  with  the 
environment,  were  set  to  a  very  high  value  of  700  W/m2  °C.  Initial  mesh  temperatures  were 
set  at  30°  C  and  the  environmental  temperature  was  maintained  at  26°  C.  After  100  time 
steps,  the  boundary  temperature  at  z=0  was  reset  to  20°  C  and  the  program  was  run  for 
100,000  time  steps,  1  second  each.  Due  to  the  high  value  of  the  heat  transfer  coefficients 
used,  a  rather  flat  temperature  profile  of  26°  C  was  established  and  maintained  throughout 
the  mesh  except  for  a  short  drop  to  20°  C  visualized  near  the  base  of  the  cylinder,  as  is  to 
be  expected. 

In  test  #6,  blood  perfusion  was  activated  at  the  basal  values  listed  in  Table  4.  Other 
parameters  were  maintained  inactive.  Initial  mesh  and  arterial  temperatures  were  set  at  20° 
C.  After  the  initial  1000  time  steps,  0.1  second  each,  both  mesh  and  arterial  temperatures 
were  reset  to  30°  C  at  z=0.  The  test  was  run  for  a  total  of  200,000  time  steps  and  mesh 
temperatures  converged  on  30°  C  and  remained  stable  for  the  duration  of  the  test. 

A  similar  test  was  run  with  the  addition  of  heat  exchange  between  the  major  blood 
vessels  and  the  mesh  points.  Results  of  this  test  #7  were  essentially  similar  to  those  of  test 
#6. 


In  test  #8  counter-current  heat  exchange  between  the  major  blood  vessels  themselves 
was  also  activated.  Running  conditions  were  identical  to  those  of  test  #6  except  that  a  total 
of  1 ,000,000  time  steps  were  utilized.  Mesh  temperatures  have  stabilized  at  30°  C  after 
100,000  time  steps. 

As  was  to  be  expected,  the  execution  of  the  code  was  sensitive  to  the  size  of  the  time 
step  and  the  number  of  spatial  divisions  used  in  the  numerical  code.  These  two  topics  are 
discussed  separately.  It  is  firstly  noted  that  the  method  of  alternating  directions  applied  to  the 
solution  of  the  mesh  temperatures  yields  an  unconditionally  stable  scheme  of  solution  [10]. 
Thus,  the  source  for  this  sensitivity  must  reside  in  Equations  (7)  and  (8)  representing  the  heat 
balance  in  the  major  blood  vessels.  These  two  equations  are  essentially  first  order  ordinary 
differential  equations.  Thus,  the  terms  multiplying  the  independent  variables  on  the  right 
hand-sides  may  be  used  to  estimate  the  maximal  time  step  that  will  ensure  stability  of  the 
Euler's  scheme  used  to  solve  them. 

In  the  present  analysis  this  maximal  time  step  is  determined  by  calculating  the 
numerical  values  of  the  coefficients  of  the  independent  variables  in  Equations  (41 )  and  (42). 
The  larger  of  the  two  values,  TOTAL,  is  then  substituted  into  the  following  equation: 


A  x 


0.4 

TOTAL 


(50) 
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to  yield  the  time  step  which  ensures  numerical  stability.  Values  obtained  by  Equation  (50)  are 
conservative  since  a  factor  of  2  may  be  used  instead  of  0.4  [21].  Experience  with  running  of 
the  code  proved  that  this  requirement  on  the  time  step  could,  indeed,  be  relaxed  somewhat 
without  adversely  affecting  the  stability  of  the  code. 

The  sensitivity  of  the  code  to  the  number  of  divisions  used  in  the  numerical  network 
became  apparent  when  an  overall  steady  state  heat  balance  was  calculated  for  the  finger 
model.  In  all  cases  studied  the  number  of  divisions  in  the  radial  direction  was  kept  constant 
at  12  (Table  6).  A  relatively  simple  combination  of  parameters  was  used  in  the  computations 
in  which  the  finger  was  assumed  to  be  insulated  from  the  environment,  no  metabolic  heat 
was  generated  and  no  counter-current  heat  exchange  between  the  major  blood  vessels  was 
allowed.  In  addition,  the  thermal  conductivities  of  all  tissue  compartments  were  made  uniform 
at  the  value  of  the  muscle.  Under  these  conditions,  the  only  heat  supply  to  the  tissue  was 
due  to  blood  perfusion  and  the  only  heat  removal  mechanism  was  by  conduction  at  the  base 
of  the  finger,  i.e.,  at  z=0. 

Figure  5  shows  the  ratios  calculated  for  the  heat  transported  by  the  major  blood 
vessels  to  the  heat  conducted  away  as  a  function  of  the  number  of  axial  divisions.  It  is 
evident  that  a  heat  balance  is  not  satisfied  for  the  smaller  number  of  divisions  in  the  axial 
direction.  Only  when  20  divisions  are  used,  the  heat  source  essentially  equals  the  heat  sink 
to  satisfy  a  heat  balance. 

A  similar,  but  more  involved,  set  of  benchmark  tests  was  also  run.  In  this  set  both 
metabolic  heat  production  and  heat  exchange  with  the  environment  at  the  finger  tip  were 
included  in  addition  to  blood  perfusion  and  heat  conduction  at  the  finger  base.  A  steady-state 
energy  balance  offset  of  about  32%  was  initially  obtained  for  25  axial  and  9  radial  divisions. 
This  offset  gradually  dropped  to  less  than  1%  when  75  axial  divisions  were  used.  This  value 
was  deemed  quite  satisfactory  for  a  steady-state  energy  balance  and  served  as  an  additional 
verification  of  the  numerical  code. 

Temperature  distributions  along  the  external  surface  of  the  model  are  plotted  in  Fig. 
6  also  as  functions  of  the  number  of  axial  divisions.  It  is  clearly  seen  that  the  final 
temperature  obtained  is  a  function  of  the  number  of  divisions  in  the  axial  direction.  For  the 
particular  case  studied  here  it  seems  that  20  divisions  yield  quite  an  accurate  result.  This, 
however,  required  a  much  longer  running  time  than  for  the  fewer  divisions.  Thus,  the  desired 
accuracy  of  the  results  may  have  to  be  determined  by  considerations  such  as  the  total  CPU 
time  required  for  running  the  program  for  a  given  computer. 

Based  on  the  series  of  benchmark  tests  as  outlined  here,  it  may  be  concluded  that  the 
code  written  for  the  model  is  stable  and  converges  to  reasonable  results  for  the  entire  range 
of  parameters  considered  here. 
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LENGTH  OF  CYLINDER,  cm 

8.0 

DIAMETER  OF  CYLINDER,  cm 

1.5 

DIAMETER  OF  ARTERY,  cm 

DIAMETER  OF  VEIN,  cm 

DISTANCE  BETWEEN  ARTERY  AND  VEIN,  cm 

Eq.  (B9) 

HEAT  TRANSFER  COEFFICIENT  BETWEEN  ARTERY 

OR  VEIN  AND  THE  TISSUE  (Ua  or  Uv) 

Eq.  (B6) 

COUNTERCURRENT  HEAT  EXCHANGE  COEFFICIENT 
BETWEEN  ARTERY  AND  VEIN  (h,  or  hv) 

Eq.  (B8) 

NUMBER  OF  DIVISIONS  IN  THE  RADIAL  DIRECTION 

CORE 

MUSCLE 

FAT 

SKIN 

TOTAL 

3 

3 

2 

4 

12 

Table  6:  Parameters  used  in  the  numerical  computations. 
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Figure  5  (top):  Ratios  of  steady  state  heat  flow  in  vs.  heat  flow  out  as  affected  by 

the  number  of  divisions  in  the  axial  direction. 


Figure  6  (bottom):  Steady  state  temperature  distributions  of  the  external  surface  of 
the  finger  model  as  affected  by  the  number  of  divisions  in  the  axial 
direction. 


RESULTS  AND  DISCUSSION 


A  number  of  cases  are  considered  to  demonstrate  the  range  of  capabilities  of  this 
model.  Additional  parameters  used  in  these  demonstrations  are  listed  in  Table  6. 

Figures  7-10  show  the  steady  state  temperature  distributions  in  the  finger  model  for 
combinations  of  insulation  (bare  vs.  gloved  finger),  wind  velocities  and  finger  blood  flow.  In 
all  these  cases  the  environmental  temperature  was  maintained  at  0°  C.  and  finger  base 
temperature  as  well  as  incoming  arterial  blood  temperature  were  kept  constant  at  30°  C. 

All  four  figures  demonstrate  the  major  role  played  by  blood  flow  in  the  thermal 
economy  of  the  finger.  It  is  clearly  noted  that  rather  comfortable  temperatures  are  maintained 
in  the  finger  for  as  long  as  blood  flow  remains  high  (at  basal  level  in  this  case).  The 
exception  is  the  case  of  the  bare  finger  in  a  windy  environment  of  15  km/hr  in  which  the 
enhanced  heat  loss  offsets  much  of  the  beneficial  effect  of  high  blood  flow  to  the  finger. 

In  all  cases  studied,  temperature  of  the  distal  segments  of  the  finger  dropped  to  very 
low  levels  and  almost  equilibrated  with  the  environment.  This  is  the  case  even  when  a  two- 
layer  glove  is  donned  on  the  hand  as  is  also  suggested  in  another  study  [22].  The  only 
difference  among  the  cases  studied  here  is  in  the  time  course  of  change  in  these 
temperatures.  This  difference  is  shown  in  Fig.  1 1  in  which  finger  tip  skin  temperatures  are 
plotted  vs.  time  for  all  4  cases.  Low  blood  flow,  at  the  nutritional  level,  which  is  to  be 
expected  for  this  low  environmental  temperature,  was  assumed  for  all  cases. 

As  seen  in  Fig.  11,  the  bare  finger  in  windy  air  will  be  the  quickest  to  drop  in 
temperature.  It  would  practically  equilibrate  with  the  environment  after  about  10  minutes.  The 
gloved  finger,  under  the  same  windy  environment,  would  be  much  better  protected  and  would 
require  about  60  minutes  before  it  equilibrates  with  the  environment.  Interestingly,  a  bare 
finger  in  still  air  seems  to  be  better  protected  than  a  gloved  finger  in  windy  air. 

These  results  may  also  be  presented  in  terms  of  endurance  times,  defined  as  the  time 
for  any  temperature  on  the  finger  to  reach  5°  C.  Accordingly,  the  endurance  times  for  the 
bare  and  gloved  fingers  in  windy  air  would  be  about  2.5  and  22  min,  respectively.  These 
times  would  be  longer,  at  32  and  43  minutes,  respectively,  when  the  bare  and  gloved  fingers 
are  exposed  to  still  air. 
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TEMPERATURE. 


Figure  7:  Steady  state  temperature  distributions  of  a  bare  finger  in  still  air  for  basal 

and  nutritional  blood  flows. 
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TEMPERATUR 


Figure  8:  Steady  state  temperature  distributions  of  a  gloved  finger  in  still  air  for 

basal  and  nutritional  blood  flows. 
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Steady  state  temperature  distributions  of  a  bare  finger  in  windy  air  for 
basal  and  nutritional  blood  flows. 


TEMPERATURE, 


Figure  10:  Steady  state  temperature  distributions  of  a  gloved  finger  in  windy  air  for 
basal  and  nutritional  blood  flows. 
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TEMPERATURE,  °C 


60 


Figure  1 1 :  Temperature  variations  on  the  dorsal  tip  of  the  bare  and  gloved  finger 
model  with  nutritional  blood  flow  for  still  and  windy  air. 


-  32  - 


The  effect  of  counter-current  heat  exchange  between  the  major  blood  vessels  is 
demonstrated  in  Fig.  12.  Numerical  values  of  the  parameters  used  for  this  figure  are  listed 
in  Table  6.  The  two  groups  of  curves  in  Fig.  12  show  the  steady  state  arterial  and  venous 
temperature  distributions  along  the  finger  with  and  without  counter-current  heat  exchange 
between  these  vessels.  As  soon  as  this  mechanism  is  activated,  the  arterial  temperature 
seems  to  drop  considerably  due  to  the  exchange  of  heat  with  the  cooler  vein. 

The  main  purpose  of  counter-current  heat  exchange  is  to  conserve  body  heat  in  a  cold 
environment.  This  is  effected  by  firstly  depriving  the  extremity  of  the  rich  supply  of  blood,  as 
is  assumed  here  by  dropping  blood  flow  from  basal  to  nutritional  level.  An  additional  effect 
is  achieved  by  lowering  the  temperature  of  the  extremity  through  the  thermal  coupling  which 
exists  between  the  major  blood  vessels  and  the  tissues.  The  arterial  temperature  along  the 
extremity  is  made  to  loose  heat  by  counter-current  heat  exchange  to  the  cooler  vein.  This, 
in  turn,  causes  tissue  temperatures  to  decrease  as  the  artery  constitutes  the  main  heat 
source  for  the  extremity.  In  the  case  shown  in  Fig.  12,  about  0.093  W  is  lost  to  the 
environment  from  the  finger  which  decreases  to  0.08  W  for  a  reduction  of  about  14%  in 
finger  heat  loss  when  counter-current  heat  exchange  is  activated. 

Another  case  of  cold  induced  vasodilatation  (CIVD)  in  the  finger  is  shown  in  Fig  13. 
CIVD  is  known  to  occur  in  a  percentage  of  the  population  and  is  manifested  by  rather 
periodic  increases  in  finger  skin  surface  temperatures,  e.g.,  [23,24].  Although  the  precise 
mechanism  for  this  phenomenon  is  not  thoroughly  understood,  there  is  ample  evidence  to 
indicate  that  intermittent  increases  in  the  otherwise  constricted  blood  flow  to  the  finger  cause 
these  temperature  changes. 

In  calculating  the  data  for  Fig.  1 3,  it  was  assumed  that  the  periodic  blood  flow  changes 
may  best  be  approximated  by  triangular-shaped  surges.  These  surges  were  assumed  to 
occur  only  in  and  adjacent  to  the  tip  of  the  finger  while  blood  supply  to  the  other  segments 
of  the  finger  remained  unchanged.  The  initial  temperature  of  the  entire  bare  finger,  exposed 
to  still  air  at  0°  C,  was  set  at  30°  C.  At  the  beginning  of  the  exposure,  blood  flow  in  the  finger 
was  assumed  to  drop  to  the  nutritional  level  (Table  4).  This  situation  was  assumed  to  persist 
for  20  min.  Next,  a  3  minute  linear  increase  to  10  times  the  initial  value  in  blood  flow  to  the 
tip  of  the  finger  was  allowed  followed  by  a  symmetrical  decrease  back  to  the  nutritional  level. 
Nutritional  blood  flow  was  next  maintained  for  15  minutes  and  was  followed  again  by  an 
identical  triangular-shaped  change  in  blood  flow  to  the  tip  of  the  finger.  During  the  final  15 
minutes  of  the  run,  blood  flow  was  reset  to  the  nutritional  level. 

Skin  temperature  variations  are  shown  in  Fig.  13  for  three  locations  along  the  finger. 
The  curves  are  plotted  for  the  base,  the  middle  point  along  the  finger  and  for  the  tip  of  the 
finger.  The  solid  line  represents  finger  tip  temperature  variations  for  constant  nutritional  blood 
flow  without  the  periodic  bouts  of  CIVD.  It  thus  provides  a  worst  case  scenario  for 
comparison  purposes.  It  is  clearly  seen  that  increases  in  blood  flow  causeincreases  in  finger 
temperatures,  as  is  to  be  expected.  These  changes  are  more  pronounced  at  the  tip  of  the 
finger  than  at  the  more  proximal  locations  primarily  because  blood  flow  changes  due  to  CIVD 
are  assumed  to  take  place  in  this  area  only. 
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Figure  1 3:  Dorsal  temperature  variations  at  the  tip  and  middle  of  the  finger  model  for 

cold  induced  vasodilatation. 


-  35  * 


BLOOD  PERFUSION,  kg/m 


Another  interesting  result  relates  to  the  course  of  change  in  finger  tip  temperature  for 
the  case  shown  here.  It  seems  that  CIVD,  which  may  be  characterized  as  a  heating  source, 
causes  the  tissue  temperatures  to  increase  noticeably  for  as  long  as  it  is  active.  Once  this 
mechanism  gets  shut  off,  tissue  temperatures  resume  the  exponential-like  decay  to  levels 
close  to  those  attained  without  CIVD.  This  decay  is  enhanced  by  the  larger  temperature 
difference  between  the  tissue  and  the  environment  that  has  been  established  as  a  result  of 
CIVD.  As  a  matter  of  fact  the  temperature  difference  at  the  finger  tip  after  one  hour  between 
the  case  involving  CIVD  and  the  one  without  it  would  be  a  mere  0.8°  C,  for  the  case 
presented  here. 

It  is  recognized  that  the  details  shown  in  Figs.  7-13  are  dependent  on  the  parameters 
and  assumptions  used  in  the  computations.  However,  certain  trends  are  clearly  indicated 
which  will  likely  vary  in  detail  and  magnitude  as  the  values  of  these  parameters  are  altered. 
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GLOSSARY 


VARIABLES 

a  -  dimensionless  ratio  of  cylinder  length  to  radius  multiplied  by  the  ratio  of  blood  to 
tissue  thermal  diffusivities,  Eq.  (23). 

A(.ofv)  *  dimensionless  ratio  of  cylinder  to  blood  vessel  radii,  Eq.  (43). 

[AJ  -  coefficient  matrix  in  radial  direction,  Eqs.  (37)  and  (38). 

[AJ  -  coefficient  matrix  in  axial  direction,  Eqs.  (37)  and  (38). 

B(a  or  v)  ‘  dimensionless  ratio  of  blood  flow  rate  to  blood  mass  contained  in  a  vessel 
element  per  unit  of  normalized  time,  Eq.  (44). 

Bi  -  Biot  modulus  indicating  the  dimensionless  ratio  between  heat  convected  by  the 
environment  to  heat  conducted  in  the  cylinder,  Eqs.  (28)  and  (29). 

c  -  specific  heat,  J'kg'VC1  . 

h  -  heat  transfer  coefficient  between  the  cylinder  and  the  environment  at  the 
circumferential  surface,  W*m'2*°C 1 

h1  -  heat  transfer  coefficient  between  the  cylinder  and  the  environment  at  the  cylinder 
tip,  W*m  2«°C 1 

h(r  or  z)  '  dimensionless  distance  between  two  adjacent  nodes  in  the  radial  or  axial  direction, 
respectively. 

hav  -  heat  transfer  coefficient  between  concomitant  artery  -  vein  pairs,  W^C1  . 

H(«v  or  va>  -  dimensionless  heat  transfer  coefficient  between  concomitant  artery  -  vein  pairs,  Eq. 
(46). 

k  -  thermal  conductivity,  W«m'1«°C'1 

L  -  cylinder  length,  m. 

m,aorv)  -  mass  flow  rate  of  arterial  or  venous  blood,  kg*s\ 

M  -  total  number  of  nodal  points  in  axial  direction. 

M(aorv)  ‘  mass  of  arterial  or  venous  blood  contained  in  a  vessel  element,  kg. 
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N  -  total  number  of  nodal  points  in  radial  direction. 
qm  -  volumetric  metabolic  heat  generation  rate,  W«m 3 . 
q  -  dimensionless  volumetric  heat  generation  rate,  Eq.  (25). 
r'  -  radial  coordinate,  m. 

r  -  dimensionless  radial  coordinate,  Eq.  (19). 

R  -  cylinder  radius,  m. 

{S}  -  vector  of  terms  in  matrix  equations,  Eqs.  (37)  and  (38). 

SUM2  l  '  numerical  summations,  Eqs.  (47)  -  (49),  respectively. 

SUM3 ) 

t*  -  time,  s. 

X  -  temperature,  °C. 

Tj  -  temperature  at  base  node,  °C. 

T*  -  initial  temperature  distribution  in  the  cylinder,  °C. 

Temp  -  reference  temperature,  °C. 

T  -  dimensionless  temperature,  Eq.  (22). 

u(aorv>  •  heat  transfer  coefficient  between  an  artery  or  a  vein,  respectively,  and  the 
surrounding  tissue,  W-m'^C'1  . 

U(aorv)  •  modified  dimensionless  heat  transfer  coefficient  between  an  artery  or  a  vein, 
respectively,  and  the  surrounding  tissue,  Eq.  (45). 

wb  -  volumetric  blood  perfusion  rate,  kg-m'^s'1 . 

W  -  dimensionless  volumetric  blood  perfusion  rate,  Eq.  (26). 

Y(aorV)  -  dimensionless  heat  transfer  coefficient  between  an  artery  or  a  vein  and  the 
surrounding  tissue,  Eq.  (27). 

z  -  axial  coordinate,  m. 

z  -  dimensionless  axial  coordinate,  Eq.  (20). 
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GREEK  LETTERS 


a  •  thermal  diffusivity,  m2*s'\ 

Y  -  ratio  of  radial  nodal  divisions  immediately  preceding  to  immediately  following  a  tissue 
compartment  interface,  Eq.  (A.10). 

5  -  dimensionless  half  a  time  step,  Eq.  (34). 

p  -  density,  kg»m'3 . 

t  -  dimensionless  time,  Eq.  (21). 

4*  -  dimensionless  ratio  of  blood  to  tissue  thermal  inertias,  Eq.  (24). 

SUPERSCRIPTS 

*  -  dimensional  quantity. 

-  -  average  value. 

SUBSCRIPTS 

a  -  arterial, 

b  -  blood. 

i  -  integer  (radial  direction), 
i  -  initial. 

j  -  integer  (axial  direction), 

n  -  integer  (time). 


v  -  venous. 


z  -  axial. 

0  -  environmental. 
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APPENDIX  A  -  Derivation  of  the  equations  for  the  various  matrix  nodal  points 


In  this  section  the  elements  in  the  matrix  equations,  (37)  and  (38),  are  derived.  Derivation 
is  performed  for  the  different  nodal  points  included  in  the  domain  for  which  the  solution  is 
sought,  i.e.,  center  nodes,  interface  nodes,  etc.  shown  in  Fig.  A.1 .  The  derivation  begins  with 
rewriting  the  general  discretized  partial  differential  equation,  Equation  (35): 
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where  Tav  is  a  time-average  tissue  temperature  given  generally  by: 


(A.  2) 


The  specific  conditions  at  the  various  nodal  points  are  now  substituted  into  Equation  (A.l ). 


Center  i-node;  regular  j-node 

In  order  to  satisfy  the  adiabatic  condition  specified  for  this  boundary,  Equation  (14),  and 
to  retain  a  truncation  error  of  the  order  of  0(h*),  a  central-difference  numerical  approximation 
is  employed: 
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Figure  Ai:  Identification  of  nodal  points  in  the  numerical  grid 
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with  geometrical  details  shown  in  Fig.  A.2.  It  follows  from  Equation  (A.3): 

Ti+1  =  =  r0  =  0  (A. 4) 


which  yields  upon  substitution  into  Equation  (A.1): 
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External  i-node;  regular  j-node 

The  boundary  condition  to  be  satisfied  at  these  nodes  is  Equation  (15).  The  numerical 
approximation  to  this  equation,  depicted  in  Fig.  A.3,  which  retains  the  desired  truncation 
error,  is: 

ff  “  =  B1  (T°  "  Tl)  (A,6) 


which  yields: 
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Next,  a  correction  is  applied  to  the  equivalent  tissue  temperature  at  the  external  node  to 
account  for  the  difference  in  the  heat  exchange  with  capillary  perfusion  [25]: 
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radial  divisions 


Figure  A2:  Schematic  diagram  of  the  center  node 


r  =  1 


Figure  A3:  Schematic  diagram  of  the  external-regular  node 
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Substitution  of  Equations  (A.7)  and  (A.8)  into  Equation  (A.1)  obtains: 
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Interface  i-node;  regular  j-node 

(including  a  regular  i-node;  regular  j-node  by  setting  y=1) 


Interface  nodes  define  the  boundaries  between  the  various  compartments,  or  organs,  of 
the  limb  model,  e.g.,  core-muscle  interface,  etc.  Numerically  the  subdivisions  in  each  of  the 
compartments  may  be  different  depending  on  the  extent  of  any  specific  compartment  and  the 
desired  numerical  details.  Additionally,  tissue  properties  generally  vary  among  compartments 
which  further  justifies  the  definition  of  these  interfaces.  Figure  A.4  depicts  the  interface  region 
between  two  adjacent  tissue  compartments  in  the  radial  direction.  The  ratio  of  the  sub¬ 
division  immediately  preceding  the  interface  to  that  following  it,  is  defined  by: 

Y  =  — ^  (A.  10) 

Ajt 

for  simplicity,  the  superscripts  in  Equation  (A.10)  are  omitted  with  the  convention: 
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Figure  A4:  Schematic  diagram  of  the  interface  node 
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The  following  numerical  approximations  are  used  in  order  to  retain  a  truncation  error  of  order 
O(hf)  [10]: 
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After  Equations  (A.10)  -  (A.13)  are  substituted  into  Equation  (A.1),  the  following  expression 
is  obtained: 
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End  j-node;  regular  i-node 

In  the  present  analysis  it  is  assumed  that  the  major  blood  vessels  terminate  (artery)  and 
originate  (vein)  in  the  j*node  immediately  preceding  the  end  j-node,  as  depicted  in  Fig.  A.5. 
Thus,  heat  exchange  with  the  circulatory  system  is  performed  at  this  node  with  the  capillary 
bed  only.  The  boundary  condition  at  this  node  is  Equation  (17)  which  is  approximated  by: 

f  ~  ^2*,—  -  “»<*«  -V  (A-15) 


-  A8- 


CAPILLARY 

PERFUSION 


jre  A5:  Schematic  diagram  of  the  end  j  -  regular  i  node 


,9- 


which  is  rewritten  as: 
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End  j-node;  external  i-node 

At  this  node,  which  is  depicted  in  Fig.  A.6,  Equations  (15),  (17)  {or  their  approximations 
Equations  (A.7)  and  (A.16)}  and  the  modification  given  by  Equation  (A.8)  are  to  be  satisfied. 
Additionally  the  major  blood  vessels  origination-termination  assumption  is  also  applied  to 
yield: 
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End  j-node;  center  i-node 


At  this  node  both  Equations  (A.4)  and  (A.16)  apply  to  yield: 
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End  j-node;  interface  i-node 


At  this  node  Equations  (A.  12),  (A.  13)  and  (A.  16)  are  substituted  to  yield: 
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APPENDIX  B  -  Heat  transfer  coefficients  for  the  major  blood  vessels 


Under  the  assumptions  of  the  present  study  the  two  major  blood  vessels,  an  artery  and 
a  vein,  traverse  each  body  element  in  the  axial  direction,  Fig.  1.  These  vessels  exchange 
heat  with  the  surrounding  tissue.  Additionally  under  certain  circumstances  there  may  also  be 
direct  heat  exchange  between  the  two  vessels.  Mitchell  and  Myers  [26]  assumed  that  each 
one  of  these  major  vessels  resides  alone  in  the  tissue.  An  improved  assumption  is  due  to 
Keller  and  Seiler  [27]  and  Arkin  and  Shitzer  [8]  according  to  which  an  "influence  volume"  is 
defined  by  enclosing  each  vessel  by  a  larger  concentric  cylinder,  Fig.  B.l.  The  conduction 
shape  factor  for  this  situation  is  given  by  [28] 
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where  Rc  is  the  radius  of  the  enclosing  cylinder.  As  a  first  approximation  this  radius  may  be 
set  as  the  half  distance  between  the  centers  of  the  major  blood  vessels,  D.  Thus, 
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The  amount  of  heat  exchanged  between  any  of  these  vessels  and  the  surrounding  tissue 
is  now  given  by 
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where  k  is  a  volume  weighted  average  value  for  tissue  thermal  conductivity  given  by 
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Figure  B1 :  Schematic  representation  of  an  "influence"  volume  enclosing 
a  major  blood  vessel 
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The  integral  in  Equation  (B.4)  is  performed  numerically  by  the  trapezoidal  rule.  The  overall 
temperature  difference  in  Equation  (B.5)  includes  an  average  value,  T ,  representing  the 
tissue  surrounding  the  blood  vessel.  This  average  tissue  temperature  is  calculated  by 
performing  the  integrations  indicated  in  Equations  (7)  and  (8).  For  these  integrals  a 
volummeiric  average  heat  transfer  coefficient,  ua  or  v  is  required.  This  is  obtained  by 
"distributing"  the  overall  heat  transfer  coefficient  over  the  volume  of  the  entire  body  element 
to  obtain 
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The  direct  heat  exchange  between  the  artery  and  the  vein  may  be  approximated  by 
assuming  a  conduction  shape  factor  for  two  parallel  cylinders  exchanging  heat,  the  centers 
of  which  are  spaced  by  D  [28],  Fig.  B2: 
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The  overall  heat  transfer  coefficient  between  these  cylinders  is  given  by 
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The  distance  between  the  centers  of  the  two  major  blood  vessels  is  usually  known,  or  can 
be  inferred,  for  any  body  element.  For  the  purposes  of  the  present  computations  we  assume 
this  distance  to  be  given  by  four  times  the  geometric  average  of  the  radii  of  these  vessels 
[8] 
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which  yields  upon  substitution  into  Equation  (B.8) 
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APPENDIX  C  -  Program  source  code  listing  and  operating  instructions 


Program  Environment 


1 .  Hardware  requirements  :  386  or  486  PC  with  a  80X87  math  coprocessor. 

2.  TURBO  PASCAL  6.0  installed  and  operating. 

3.  BGI  directory  loaded  in  the  TP  directory  enabling  the  graphics  interface  to  operate. 

4.  Option  Switches  - 

A.  COMPILE  menu  bar  -  DESTINATION  set  to  disk. 

B.  OPTION  menu  bar  -  COMPILER  set  to  the  8087  processing  mode. 

MEMORY  SIZE  •  STACK  SIZE  set  to  65520. 


The  program  consists  of  four  files: 

1 .  FINGER93.PAS  is  the  main  body  of  the  program. 

2.  VAR_SP.PAS  is  a  Turbo  Pascal  UNIT  containing  all  global  declarations. 

3.  TISUE_SP.PAS  is  a  Turbo  Pascal  UNIT  containing  the  procedures  necessary  to 
generate  the  coefficients  for  the  tissue  equations. 

4.  PHYS_FIN.PAS  is  a  Turbo  Pascal  UNIT  containing  all  procedures  dealing  with 
anatomical  data. 

The  system  of  Pascal  programs  utilizing  a  main  program  which  accesses  three  UNITS 
was  necessary  since  files  in  TURBO  PASCAL  6.0  cannot  exceed  65K  bytes. 
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N.B.  -  all  units  must  be  re-compiled  to  disk  each  time 
changes  are  made  to  the  code  under  consideration. 
Turbo  Pascal  requires  that  a  *.TPU  file  exist  on 
disk  for  each  UNIT  in  the  program.  Consequently, 
the  units  must  be  recompiled  each  time  the  user 
makes  changes  to  these  files.  Compiling  these  files 
with  the  destination  switch  set  to  memory  will  not 
create  the  necessary  *.TPU  files  for  the  program  to 
execute  successfully.  THE  UNITS  MUST  BE  CCMPIT.fd 


Description  of  program  Units 


UNIT  FINGER93.PAS  {  main  program  } 


1 .  Procedure  initialize; 

Generates  formula  calculations  for  all  values  required  by  the  program.  Data  for 
this  procedure  is  acquired  from  the  CONST  section  of  UNIT  VAR_SP. 

2.  Procedure  Data_Dump; 

Generates  file  output  of  initial  and  calculated  data  that  wail  be  used  by  the 
program. 


UNIT  VAR_SP.PAS 

This  unit  contains  all  CONSTant  declarations,  data  TYPEs  and  VARiable  declarations 
used  by  the  program.  For  testing  purposes,  the  initial  values  assigned  in  the  CONSTant 
section  may  be  altered.  This  area  has  been  identified  by  the  internal  documentation.  No 
changes  should  be  made  in  the  TYPE  sections  or  the  VAR  sections  since  these  will  have  a 
substantial  impact  on  program  execution.  VAR_SP  must  be  saved  and  recompiled  to  disk 
each  time  the  test  values  are  changed.  This  unit  contains  six  sub-modules  which  function 
as  utilities  of  calculation  procedures  for  the  generation  of  coefficients  required  by  the  tissue 
equations. 
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1.  FUNCTION  GAMMA 

Calculates  the  gamma  values  as  required. 

2.  PROCEDURE  PAUSE_CRT 

A  utility  that  allows  the  user  to  halt  the  programs  execution  as  required. 

3.  PROCEDURE  DERIVE_CENTER_TEMP 

Derives  the  temperature  at  I  =  0. 

4  PROCEDURE  DUMP_HALF_ARRAYS 

Outputs  the  all  values  of  the  temperature  arrays  at  the  half  time  step. 

5.  PROCEDURE  DUMP_T_ARRAYS 

Outputs  the  all  values  of  the  temperature  arrays. 

6.  PROCEDURE  CALC_TA_TV_HALF 

Calculates  the  values  for  arterial  and  venous  temperatures.  Variables  TAJ3AR, 
TVJ3AR,  TA_NODE,  and  TV_NODE  are  calculated  here  in  preparation  for  the 
calculation  of  the  S  VECTOR  values. 


UNIT  TISUE_SP.PAS 

1 .  Procedure  Axial_Dir; 

Imports  the  I  and  J  coordinates,  identifies  the  type  of  node  using  enumerated 
data  types  and  calculates  three  coefficients  for  the  given  l,J  location.  These  are 
stored  in  the  ORIG_A_Z_COEF  array  which  is  used  for  tissue  traversal  in  the 
Radial  Direction  to  generate  the  constant  term  prior  to  accessing  the  code  for 
Thomas’  Algorithm. 

2.  Procedure  Rad_Dir; 

Generates  coefficients  utilizing  the  same  approach  described  above.  The 
coefficients  are  stored  in  array  ORIG_A_R_COEF  and  are  used  in  the  axial 
traversal  of  the  tissue  matrix  to  generate  the  constant  term  prior  to  accessing 
the  code  for  Thomas'  Algorithm. 


3.  Procedure  Establish_coef_arrays; 

Driver  block  containing  nested  loops  which  call  procedures  1  and  2  listed 
above,  thereby  loading  the  coefficient  arrays. 
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4.  Procedure  Traverse; 

Contains  the  following  nested  procedures: 

a.  Procedure  Thomas; 

Contains  code  for  Thomas  Algorithm. 

b.  Procedure  Flrst_Traverse_Radial; 

Moves  appropriate  coefficients  from  arrays  generated  earlier  into  variables 
needed  by  Procedure  Thomas.  This  procedure  then  calls  the  Thomas 
algorithm  procedure  and  places  the  new  temperatures  into  the  tissue 
temperature  matrix.  A  single  dimension  array  is  used  as  an  intermediary 
data  structure  in  order  to  hold  the  l,J  coordinates  and  corresponding 
temperatures. 

c.  Procedure  Second_Traverse_Axial; 

Organized  in  the  same  manner  as  the  First_Traverse  procedure.  This 
procedure  traverses  the  temperature  matrix  in  the  Axial  direction. 


d.  Main  Block  of  Procedure  Traverse; 

Driver  block  for  establishing  the  temperature  matrix  and  nested  looping 
structure  for  the  number  of  time  steps  and  print  intervals. 


UNIT  PHYS  SP.PAS 


This  unit  contains  the  following  five  procedures: 


1 .  PROCEDURE  Anatomical Jnput; 

2.  PROCEDURE  Num_Di visions; 

These  procedures  load  variables  with  the  necessary  anatomical  data. 

3.  PROCEDURE  Width.Calculation; 

Determines  the  width  of  each  organ,  H-  values  and  gamma  values. 


4.  PROCEDURE  Estjnterfaces; 

Establishes  the  location  of  the  organ  interfaces  and  stores  them  in  a  data 
structure  of  type  SET. 


5.  PROCEDURE  Create_Arrays; 

Establishes  a  set  of  values  for  each  organ  and  a  series  of  arrays  containing  R_i 
and  Hr-  values  indexed  to  each  I  location. 
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OUTPUT  FILES 


The  following  two  external  output  files  are  automatically  created  when  the  program  is 
executed. 

t.  DVALUES.DAT 

File  contains  a  listing  of  all  significant  values  used  during  the  current  program 
run.  Some  of  these  are  echo  prints  of  values  entered  in  the  CONST  section  of 
VAR_SP.PAS,  while  others  are  the  result  of  formula  calculations. 

2.  R_TEMPS.DAT 

This  file  contains  a  printout  of  the  time  step,  arterial,  venous,  and  tissue 
temperatures  generated  of  at  each  time  step,  at  user  determined  intervals 
during  program  execution. 


OVERVIEW  OF  PROGRAM  OPERATION 


1 .  Prior  to  running  the  program,  be  certain  that  the  following  programming  conditions 
have  been  established  : 

A.  Turbo  Pascal  has  been  properly  installed  with  all  the  BGI  files  in  the  BGI 
directory  within  the  TP  directory. 

B.  COMPILE  menu  bar  DESTINATION  has  been  set  to  DISK. 

C.  OPTIONS  menu  bar  has  been  u  o  set  the  compiler  to  the  8087  numeric 

processor. 


2.  Changing  program  parameters: 

A.  Accessible  Parameters  are  found  in  the  VAR_SP.PAS  program  unit.  These 
have  been  separated  from  the  rest  of  the  code  by  the  following  sets  of 
comments: 
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(*****************************************************************) 
(*******************  g^D  CGNIHQL  POINT  SECTION  *******************) 
(***********  NOTHING  SHOULD  BE  ALTERED  BELOW  THIS  POINT  **********) 
(***********★*****************************************************) 


B.  The  parameters  within  the  CONTROL  POINT  SECTION  include: 
GC_DELTA_SEC  {  delta  value  in  seconds  } 

GCJ.ENGTH  and  GC_DIAMETER  { in  centimeters  } 
Radial  Divisions 
Axial  Divisions 

TA,  TV,  TISSUE_TEMP,  ENVIRONMENTAL_TEMP,  and 
NORMALIZING_TEMP  { in  Degrees  C.} 

BASAL  METABOLIC  RATES 

BASAL  BLOOD  FLOW  RATES 


VARIABLES  IN  OTHER  SECTIONS  OF 
THE  CODE  MUST  NOT  BE  ALTERED! 


3.  When  changes  are  made  to  any  of  the  above  listed  parameters,  the  program  must 
be  re-compiled  to  disk  before  execution  in  order  to  generate  the  appropriate  *.TPU 
files. 
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N.B.  -  Turbo  Pascal  imposes  a  65K  constraint  with  respect  to  variable  memory 
utilization.  A  compiler  error  will  occur  if  the  radial  and  axial  division 
entries  generate  an  excessively  large  tissue  matrix. 

{  ERROR  96  -  TOO  MANY  VARIABLES  } 

To  avoid  this  problem  it  is  recommended  that  the  following 
limitations  be  observed  : 

TOTAL  AXIAL  DIVISIONS  <=  12 

TOTAL  RADIAL  DIVISIONS  <=  10 


4.  DELTA  VALUE  -  The  access  point  for  delta  is  through  variable  GC_DELTA_SEC 
which  is  usually  represented  as  a  fraction  of  a  second.  The  program  uses 
GC_DELTA_SEC  to  calculate  a  USER  DETERMINED  normalized  delta  value  for 
program  execution.  The  code  also  calculates  a  PROGRAM  DETERMINED 
normalized  delta  value  which  is  the  MAXIMUM  value  for  delta  based  on  the  existing 
program  parameters.  The  USER  DEFINED  value  must  be  less  than  or  equal  to  the 
PROGRAM  DETERMINED  maximum  delta  value  for  the  program  to  function 
effectively.  However,  if  time  constraints  are  important,  the  USER  DEFINED  delta 
value  may  slightly  exceed  the  program  determined  maximum  delta.  The  numerical 
integrity  of  the  program  results  are  directly  influenced  by  the  size  of  the  USER 
DEFINED  delta  value.  There  will  be  an  obvious  breakdown  in  the  numeric  output 
if  the  USER  DEFINED  delta  value  is  too  large. 


5.  APPROXIMATE  RUNNING  TIME  : 


Hardware: 

Axial  Divisions: 
Radial  Divisions: 
Running  Time: 


486  PC  running  at  33  MHz 
10 
13 

20,000  time  steps  takes  approximately  18  minutes. 


-  C7  - 


Running  the  Program 


1 .  Be  sure  that  the  TURBO  PASCAL  6.0  environment  has  been  properly 

configured. 

a.  Compiler  and  Option  switches  have  been  properly  set.  (see  Program 
Environment.) 

b.  FINGER93.PAS,  VAR_SP.PAS,  PHYS_SP.PAS,  TISUE.PAS  files  currently 
reside  in  the  TP  directory. 

c.  FINGER93.PAS,  VAR_SP.PAS,  PHYS_SP.PAS,  TISUE.PAS  files  have  been 
opened  on  the  TURBO  PASCAL  DESKTOP. 

2.  If  any  of  the  CONTROL  POINT  parameters  are  to  be  changed,  use  the  mouse  to 

double  click  on  the  VAR_SP.PAS  to  bring  the  program  text  into  the  editor. 

a.  Change  values  in  the  CONTROL  POINT  SECTION. 

b.  Compile  VAR_SP.PAS  to  disk. 

c.  If  no  compilation  errors  exist,  double  click  on  VAR_SP.PAS  to  exit  and  return 
to  the  desktop. 

3.  Single  click  on  FINGER93.PAS. 

4.  Click  RUN  on  the  menu  bar  and  click  RUN  on  the  pull  down  menu  OR  press 

<CTRL>  F9  to  execute  the  program. 

5.  Press  <RETURN>  at  the  title  screen. 

6.  Check  the  USER  DETERMINED  DELTA  value  against  the  PROGRAM 

DETERMINED  MAXIMUM  DELTA  value. 

a.  The  value  for  GC_DELTA_SEC  which  was  entered  in  the  CONTROL  POINT 
section  of  UNIT  VAR_SP.PAS  is  used  to  calculate  the  delta  value  which  will  be 
accessed  by  the  program.  The  program  also  determines  a  maximum 
acceptable  delta  value  under  the  existing  parameters.  It  is  essential  that  the 
USER  DETERMINED  DELTA  value  does  not  exceed  the  PROGRAM 
DETERMINED  MAXIMUM  DELTA  value. 

b.  At  the  time  step  prompt  enter  a  number  1. 

c.  At  the  interval  prompt  enter  a  number  1 . 
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d.  USER  DETERMINED  DELTA  and  the  PROGRAM  DETERMINED  MAXIMUM 
DELTA  values  will  appear  on  the  screen.  Be  sure  that  the  USER  DETERMINED 
DELTA  is  less  than  or  equal  to  the  PROGRAM  DETERMINED  MAXIMUM 
DELTA  value. 

e.  After  checking  the  relationship  between  the  USER  DETERMINED  DELTA  and 
the  PROGRAM  DETERMINED  DELTA  follow  the  prompts  at  the  bottom  of  the 
screen  pressing  <RETURN>  until  program  execution  terminates. 

f.  If  the  USER  DETERMINED  DELTA  is  greater  than  the  PROGRAM 
DETERMINED  DELTA,  repeat  steps  2  and  3  using  a  smaller  GC_DELTA_SEC 
value  until  the  correct  relationship  is  obtain  between  the  USER  DETERMINED 
DELTA  and  the  PROGRAM  DETERMINED  MAXIMUM  DELTA  values.  DO 

NOT  proceed  to  STEP  7  unless  a  valid  DELTA  value  has  been  obtained. 


6.  Click  on  program  FINGER93.PAS  and  then  press  <CTRL>  F9  to  begin  executing 
the  program  for  the  desired  time  interval. 

a.  At  the  time  step  prompt,  enter  the  desired  number  of  total  times  steps  for  the 
program  run.  TOTAL  TIME  STEPS  =  total  clock  minutes  *  60  / 
GC_DELTA_SEC 

b.  At  the  interval  prompt,  enter  the  time  step  interval  at  which  to  display  data  on 
the  screen  and  send  data  to  output  file  R_TEMPS.DAT 

INTERVAL  =  interval  minutes  *  60  /  GC_DELTA_SEC 


c. 


GC_DELTA_SEC  = 
total  clock  minutes  = 
print  interval  = 

TOTAL  TIME  STEPS  = 
INTERVAL  = 


0.08 

62  minutes 
1  minute 
46500 
750 


7.  When  the  program  begins  execution,  a  time  step  counter  appears  on  the  screen. 
Temperature  output  to  both  screen  and  file  R_TEMPS.DAT  will  occur  at  the  time 
step  INTERVAL  established  in  step  6.  Follow  the  prompts  at  the  bottom  of  the 
screen  to  continue  execution  or  to  view  the  graph  of  the  important  temperatures 
generated  at  the  established  INTERVAL. 
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8.  Upon  program  completion,  two  files  will  have  been  created  in  the  TP  directory: 


a.  File  D_VALUES.DAT  contains  the  CONTROL  POINT  values  used  in  program 
run  as  well  as  formula  calculations. 

b.  File  R_TEMPS.DAT  contains  the  information  that  was  printed  to  the  screen  at 
the  established  INTERVALS  during  the  program  run. 


c.  To  view  these  files  on  the  screen  use  the  MS-DOS  TYPE  command. 


PROGRAM  SOURCE  CODE  LISTINGS 


1.  FINGER93.PAS 

{ . 


program  description  :  A  numerical  solution  for  the  thermal 

exchange  of  extremities  in  a  cold  environment 
based  on  the  model  developed  by  Dr.  Avraham  Shitzer. 

written  for  :  USARIEM,  Natick 

written  by  :  Paul  Vital 
language  :  PASCAL  6.0 

WORKING  COPY  ---  FINGER  v3.0  -  AUGUST  27,  1993 


PROGRAM  FINGER_93  (INPUT, OUTPUT)  ; 

USES 

var_sp,  physjsp,  tisue_sp,  DOS,  CRT,  PRINTER; 


PROCEDURE  DESTINATION; 

{ .  determines  output  channel  printer/crt . } 

BEGIN 

CLRSCR; 

ASSIGNCRT  (F); 

REWRITE  (F); 

ASSI(2J(F2, 'D_VALUES.DAT' ) ; 

REWRITE  (F2)  ; 

ASSIGN  (OUT  DATA,  'R_TEMPS.DAT'); 

REWRITE  (CUT  DATA)  ; 

END; 


PROCEDURE  DAIAJDUMP; 

{ . . .prints  axial  and  radial  coefficients  as  well  as  results  of 

all  formula  initializations . } 


VAR 

INDEX, I, J  :  INTEGER; 

H_T0TAL  :  REAL; 

BEGIN  (*  data  dump  *) 

CLRSCR; 

WRITE  ('DUMP  DATA  [Y,  N1  ?  '); 

IF  TJPCASE  (READKEY)  =  'Y'  THEN 
BEGIN 

WRITELN  (F2,  'DATA  DUMP'  :  47)  ; 

WRITELN  (F2)  ; 

WRITE  (F2 ,  'PHYSICAL  LENGTH  =  ’ ) ; 

WRITELN  (F2,LEN  :  10:7,  'meters'  :  10); 

WRITELN  (F2,  'DIAMETER  =  '  :20,  DIAM  :  10:7, 'meters’  :  10); 
WRITELN  (F2)  ; 

WRITELN  (F2,  'Radii  of  artery  and  vein  {  units  =  meters  }'); 
WRITELN  (F2,  'R  A  =  •  :  10, R  A: 7 : 4)  ; 

WRITELN  (F2 ,  'R  V  =  '  :10,R"V: 7 : 4)  ; 

WRITELN  (F2)  ;  “ 

WRITELN (F2)  ; 

WRITELN  (F2 ,' INITIAL  VALUES - - - '); 
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WRITEIN  (F2)  ; 

WRITELN  (F2,  'TA  *  '  :25,  TA  :3)  ; 

WRITELN  (F2,  'TV  «  ':2s,  TV  :3); 

WRITELN  (F2,  'Tissue  tenperature  =  ' :25,  mesh  TEMP : 3 ) ; 

WRITELN  (F2)  ; 

WRITETN  (F2,  'Environmental  tenperature  =  ' ,  t  0  *  nornalizing_tenp:4 :2)  ; 
WRITELN  (F2,  'W  values  {  units  =  kg.  blood  /~ cu.m  tissue  /  sec  }*); 
WRITELN  (F2,  'CSRE' :28,  Core  blood  meters  :10:3); 

WRITELN  (F2,  'MUSOJ3' :28,  ttus  blood  mete.  3  :10:3); 

WRITELN  (F2,  'EAT'  :28,  fat  blood  meters: 10:3) ; 

WRITELN  (F2,  'SKIN'  :28,  skin  blood_meters  :10:3) ; 

WRITELN  (F2)  ; 

WRITELN  (F2,  'organ  radius  to  finger  ratio  {  R  i  /  R  }'); 

WRITELN  (F2,  'CORE  RATIO' : 28,  core  ratio  :10:4); 

WRITELN  (F2,  'MUSCLE  RATIO' :28,  muscle  ratio  :10:4); 

WRITELN  (F2,  'EAT  RATIO’ : 28,  fat  ratio:10:4); 

WRITELN  (F2,  'SKUJ  RATIO' : 28,  skin  ratio: 10:4) ; 

WRITELN  (F2)  ; 

WRITELN  (F2 ,  '  thermal  conductivity  {  units  *  Watt  /  m  /  deg .  C  } ' )  ; 
WRITELN  (F2,  'CORE  THERM  0CN':28,  core  therm  con  :10:3); 

WRITELN  (F2,  'MUSCLE  THER  M  OCN' :28,  mus  therm  con  :10:3); 

WRITELN  (F2,  'EAT  THERM  CCN' :28,  fat  therm  con  :10:3); 

WRITELN  iF2,  'SKIN  THERM  OCN'  : 28,  skln_therm  con  :10:3); 

WRITELN  (F2,  'BLOOD  THERM  OCN'  28,  blood  therm  con  :10:3); 

WRITELN  (F2)  •  ™ 

WRITELN  (F2/'heat  edacity  {  units  =  J  /  kg  /  deg.  C  } ' )  ; 

WRITELN  (F2,  'CORE  HEAT  CAP' :28,  core  heat_cap  :  10:2); 

WRITE" N  (F2,  'MUSCLE  HEAT  CAP' : 28, rrus_fieat_cap  :10:2); 

WRITLN  (F2 ,  'EAT  HEAT  CAP' : 28,  fat  heat_cap:  10:2)  ; 

WRITELN  (F2,  'SKIN  HEAT  CAP'  :28,  skin  heat_cap  :10:2)  ; 

WRITELN  (F2,  'BLOOD  HEAT  CAP' : 28,  blood_heat_cap  :10:2); 

WRITELN  (F2)  ; 

WRITELN  (F2, 'density  {  units  =  kg  /  cu.  m  }'); 

WRITELN  (F2,  'CORE  DQJSITY'  :28,  core  density  :10:2); 

WRITELN  (F2,  'MCk-CLE  DENSITY'  :28,  mus  density  :10:2); 

WRITELN  (F2,  'EAT  DENSITY' :28,  fat  density  :10:2); 

WRITELN  (F2,  'SKIN  DENSITY'  : 28,  skoh  denslty  :10:2)  ; 

WRITELN  (F2,  'BLOOD  DENISTY' : 28,  blood  density  :10:2)  ; 

WRITELN  (F2)  ; 

WRITELN  (F2,  'M  =  '  ,M)  ; 

WRITELN  (F2,  'N  =  '  ,N)  ; 

WRITELN  (F2)  ; 

WRITELN  (F2,  'CORE  MUSCLE  INTERFACE  EXISTS  AT  I  VALUE  OF  ',  I_CORE_MUS)  ; 
WRITELN  (F2)  ; 

WRITELN  (F2,  'MUSCLE- FAT  INTERFACE  EXISTS  AT  I  VALUE  OF  '  ,I_MUS  EAT)  ; 
WRITELN  (F2)  ; 

WRITELN  (F2,  'FAT-SKIN  INTERFACE  EXISTS  AT  I  VALUE  OF  ' ,  I  EAT  SKIN)  ; 
WRITEIN  (F2); 

WRITELN  (F2,  'SKIN  SURFACE  EXISTS  AT  I  VALUE  OF  ' ,  SKIN  SURFACE)  ; 

WRITELN  (F2)  ; 

(*  DUMP  FORMULA  CALCULATIONS  *) 

WRITELN  (F2,  1  FORMULA  CALCULATIONS ':  50)  ; 

WRITELN  (F2)  ; 

WRITELN  (F2,  'DELTA  sec  =  '.DELTA  sec:13:10,  '  secs'); 

WRITELN  (F2,  'DELTA  =  '  .DELTA :  13  :T0,  '  normalized'); 

WRITEIN  'i 2,  'PROGRAM  CALCULATED  MAXIMUM  DELTA  VALUE  =  ',  MAX  DELTA  :  20:15) 
WRITEIN  (F2,  'USER  DETERMINED  delta  =  '  .DELTA: 20:15)  ; 

WRITELN  (F2)  ; 

WRITEIN  (F2, 'A  a  =  ',  a_a:20:10); 

WRITELN  (F2,'Av  =  ',  a_v:20:10); 

WRITELN  (F2 )  ; 

WRITE  (F2 ,  ' ******  subscript  locations  are  one  greater  than ' ) ; 

WRITELN  (F2, '  tissue  location  *****'); 

WRITEIN  (F2)  ; 

FOR  I  :=  1  TO  N  DO 
BEGIN 

WRITE  <F2,  'Q  ARRAY_PHYSC ,1, ']  =  ' ,Q  ARRAY_PHYS [I] : 15 : 10) ; 

WRITEIN  (F2,  Tq_array_nrml  [':20,i,'  T  =  ' , q_array_nrml [i] : 15:10) ; 

END; 

WRITELN  (F2)  ; 
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WRITELN  (F2,  'U  a_phys  =  1  ,U  a_phvs  :12:5,  *U  a  =  •  :10,U  a  :12:5)  ; 
WRITELN  (F2,  '0  vphys  =  ',U  vphys  -.12:5,  '0  v  =  '  :10,U  v  :12:5)  ; 

WRITELN  (F2)  ;  ~ 

VJRITELN  (F2,  1  K_AVG  =  \  K_AVG  :  20:15); 

WR  JJ  (P2)  * 

WRITEIN  (F2  /'HA  V_PHYS  =  * ,  H  A  V  PHYS  :  20:15)  ; 

WRITEIN  {F2,  'IfVv  =  h  a_v  7  12:5); 

WRITEIN  (F2,  'iTv  a  =  H  v  a  :  12:5) ; 

WRITELN  (F2)  ; 

TOR  I  :=  1  TO  N  DO 
BEGIN 

WRITE  (F2, ’w  b  nrml  [\I, •]  =  • ,w_b_nrml  [  I  )  :10:6) ; 

WRITELN  (F2,Tw_array  [‘:20,  I,  ’]  =  ',w_array  [  I  ]:10:6); 

END; 

WRITELN  (F2)  ; 

WRITELN  (F2,  ’capillary  loss  data  :’); 

WRITELN  (F2,  ’  capillary  sum  =  ' , capil  sum:20:10); 

WRITELN  (F2, '  capillary  loss  =  capillary  loss:20:10); 
WRITELN  (F2, '  starting  M_A  =  start  m_a:2$:10) ; 

WRITELN  (F2)  ; 

FOR  J  :=  1  TO  M  DO 
BEGIN 

WRITE  (F2,  'm  a  [',J, ']  =  ' ,m_atJ] :20:14) ; 

WRITELN  (F2,  Tm  v  [':20,J,  '  ]  =  ',m  v[J]:20:14); 

END; 

WRITELN  (F2)  ; 

WRITELN  (F2,  'w  b  normalizing  val  =  ' , W  b  normalizing  val  :  20:14)  ; 

WRITELN  (F2)  ; 

WRITELN  (F2,  'H  1  =  HI:  10:4); 

WRITELN  (F2,  'iTc  =  ti~C:  10:4); 

WRITELN  (F2)  ; 

WRITELN  (F2,  'B_i  C  =  ',  B_i:20:14); 

WRITELN  (F2)  ; 

FOR  INDEX  :=  1  TO  N  DO 

WRITELN  (F2,  'B  i  1  [' .index,']  =  '  ,E  I  1 [INDEX] :10:8) ; 

WRITELN  (F2)  ; 

END; 

END;  (*  data_dunp  *) 


(*=====================  FORMULAS  &  TRAVERSAL  PROCEDURES  ========; 


FUNCTION  PLACE  (  I:  INTEGER)  :  ORGAN; 

{ . translates  I  value  to  enumerated  datatype 

used  to  access  organ_values  array . } 

BEGIN 

IF  I  IN  CORE  SET  THEN 
PLACE  :=  CCORE 
ELtSE 

IF  I  IN  MUS  SET  THEN 
PLACE  :=  FfrlUSCLE 
KT.SF 

IF  I  IN  EAT  SET  THEN 

place  :=  Peat 

EL£E 

IF  I  IN  SKIN  SET  THEN 
PLACE  :=  SSKIN; 

END;  (*  place  *) 

PROCEDURE  T04P_INITTALIZE  ; 

VAR 

I,J  :  INTEGER; 

SUM  :  EXTENDED; 


BEGIN 


(*  initialize  tenperature  array  *) 

FOR  1  0  TO  N  DO 

PC®  J  :=>  1  TO  M  DO 

TEMP  GRID  II,  J]  :=  MESH  TEMP  /  NORMALIZING  TEMP; 


(*  initialize  and  normalize  venus  and  arterial  tenp  arrays  *) 

(*  assign  t  a  node  and  calculate  v  node  *) 

FOR  J  :=  1  TO  H-T  DO 

t_a_node  [J]  :*  TA  /  NORMALI Z IN3_TEMP ; 

SUM  :=  0.0; 

FOR  J  M-l  TO  M  DO 
PC®  I  :*  1  TO  N  DO 

SUM  :=  SUM  +  TEMP_GRID[I,J]  ; 

T_V_NODE [M-I]  :=  SUM  /  (N  *  2)  ; 

(*  m-2  is  original  start  point  *) 
for  j  :=  m-l  downto  1  do 

t_v_node  [j]  :=  TV  /  NORMALIZING_TEMP ; 

(*  initialize  t_v_bar,  t_a_bar  arrays  *) 

PC®  J  :=  2  TO  M-l  DO 

T  A  BAR  [J]  :=  0.5  *  (T_A  NODE  [J-l]  +  T  A  NODE  [J]  )  ; 
FOR  J  :=  2  TO  M-l  DO 

T  V  BAR  [J]  :=  0.5  *  (TV  NODE  [J-l]  +  T  V  NUDE  [J] )  ; 


(*  assign  dumny  values  to  unused  array  elements  *) 
t  a  barlm]  :=  -9999.9; 

(*  ~t“v_bar[m]  :=  -9999.9;  *) 
t  a  node[m]  :=  -9999.9; 

(*  t_v_node[m]  :  =  -9999.9;  *) 


T  A  :=  T  A  BAR; 

T^V  :=  Xy^BKR; 

old  t_a_node  :=  t_a_node; 
olcT[t_v_node  :=  t_v_node; 

END;  (*  TEMP_INITIALIZE  *) 

PROCEDURE  INITIALIZE; 

{ . initializes  all  necessary  formula  calculations 

and  tenperature  values . } 


VAR 

I,  J  ;  INTEGER; 

I  LOC  :  ORGAN; 

TOTAL,  total  a,  total_v,  SUM,  UaSQRT,  UvSQRT,  EXP1,  GAfM  :  extended; 
index  :  organ; 

function  hyperbol_cos  (x  :  extended)  :  extended; 

{  inverse  hyperbolic  cosine  function,  used  to  calculate  H_a_v  } 
var 

val,  el,e2,e3  :  extended; 

function  power  (base  :  extended;  exp  .-  integer)  .-extended; 
var 

index  :  integer; 
tenp  :  extended; 

begin 

tenp  :=  1; 
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for  index  :=  1  to  exp  do 
terrp  :>  tenp  *  base; 
power  :«  tenp; 
end;  (*  power  *) 

begin  (*  hyperbol  cos  *) 
el  :=  1/  (2  *  SQR  (x) )  ; 
e2  1  /  (4  *  POWER (x, 4) } ; 

e3  1  /  (6  *  POWER  (x,6)); 

val  :=  LN(2  *  x)  -  (1/2)  *  el  -  (3/8)  *  e2  -  (15/48)  *  e3; 
hyperbol_cos  : =  val; 
end;  (*  hyperbol_cos  *) 

BEGIN  (*  initialize  *) 


(*  INITIALIZE  ORGAN  VALUE  ARRAYS  USING  VALUES  FROM  CONST  BLOCK 

thermal  conductivity  *) 

ORG  VALS  [CCORE] .THERM  CCN  :=  core  therm_con; 

QRG~VALS  [hNUSCLE]  .THERM  CCN  :=  mus  therm  can; 

ORCrVALS  (FEAT)  .THERM  OCR  :=  fat  therm  can; 

ORG- VALS  [SSKIN]  .THERM  CCN  :=  skin  therm  con; 

ORG^VALS  [BBLOCD]  .THERHjOON  ;=  blood_therm_ccn; 

(*  heat  capacity  *) 

ORG  VALS  [CCORE]  .HEAT  CAP  :=  core  heat  cap; 

ORG- VALS  [RMUSCLE]  .HEKT  CAP  :=  mus  heat  cap; 

ORG~VALS  [FEAT]  .HEAT  CAP  :=  fat  heat  cap; 

ORG""VALS  [SSKIN]  .HEAT  CAP  :=  skin  heat  cap; 

ORG^VALS  [BBLOCD] ,HEAT_CAP  :=  blood_heat_cap; 

(*  density  *) 

ORG  VALS  [COORE]  .DENSITY  :=  core  density; 
org'Vals  [RMUSCLE]  .DENSITY  :=  mus  density; 

ORG^VALS  [FEAT] .DENSITY  :=  fat  density; 

ORG”VALS  [SSKIN]  .DENSITY  :=  skin  density; 

ORG^MLS  [BBLOCD]  .DENSITY  :=  blood_density; 

(*  basal  metabolic  rate  *) 

ORG  VALS  [COORE]  .  METAB  :=  core  me  tab; 

ORG  VALS  [NMUSCLE]  .METAB  :=  mus  me tab ; 

ORCTVALS  [FEAT]  .METAB  :=  fat  me tab; 

ORG” VALS  [SSKIN]  .METAB  :=  skin  me  tab  ; 

ORGTVALS  [BBLOCD] .METAB  :=  -9959;  (*  -9999  is  a  filler  value  *) 


(*  establish  basal  blood  flow  rate  in  cubic  me-  ers  *) 

ORG  VALS  [COORE] .METERS  BAS  BLOOD  :=  core  blood  meters ; 
ORG~VALS  [M4USCLE]  .MEIERS  BSS  BLOCD  :=  mus  blood  meters ; 
ORG- VALS  [FEAT]  .METERS  BAS  BlOCD  :=  fat  blood  meters; 
ORCTVALS  [SSKIN]  .MEIERS  BAS  BLOCD  :=  skin  blood  meters ; 
ORGTVALS  [BBLOOD]  .METERS  BAS  BLOOD  :=  -9959; 


(*  establish  A_a  and  A_v  Values  *) 

A  a  :=  SQR  (PHYS  RAD  /  R  a)  ; 

A_v  SQR  (PHYS_RAD  /  R_v) ; 

(*  establish  normalizing  value  for  w__b  *) 

W_b_normalizing_val  :=  (SQR (PHYS  RAD)  *  CM3  VALS  [BBLOOD]  .HEAT  CAP)  / 

ORG_VAES  [BBLOCD]  7THERMJXN; 

(*  establish  blood  alpha  *) 
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BL  ALPHA  :«  ORG  VALS  (BBLOQDj  THERM  CCN  / 

(ORS_VALS  [BBLOCD]  .DENSITY  *  ORG_VALS  [BBLOOD]  .  HEAT_CAP)  ; 

(*  establish  w_b_nrml_array  -  basal  blood  flow  rate  in  cubic  meters 

qarray  nrml,  b  i  1  array,  THERM  CCN  ARRAY,  HEAT  CAP  ARRAY, 
density"- ARRAY  *T  “  ~  “ 


(*  N.B.  - . - . . . - . . 

subscripts  for  tissue  data  access  represent  the  interface 
inmediately  above  the  tissue 

- . - . - . - . *) 


PC®  I  :=  1  TO  N  DO 
BEGIN 

I  DOC  :=  PLACE (I) ; 

W  ARRAY  [I]  :=  ORG  VALS  [I  ICC]  .MEIERS  BAS  BLOCD; 

Q~ ARRAY  PHYS  [I]  :=  ORG  VAES  [I  LOC]  .MEflABT 
THERM  caj  ARRAY  [I]  :=  5RG  VALSH  LOC]  .THETOI  CON; 
heat  Cap  Array  [i]  :=  orCA^ls [i~icc]  .heat  Cap,- 
DENSITY  Array  [I]  :=  ORG  ^ALS  [I~IOC]  .DENSITY; 

WbnrmT  [I]  :=  W_ARRAY[U  *  W_Bj)ORMALIZING_VAL; 

END; 

W  ARRAY  [N]  :=  0.0; 

W~b  nrml [N]  : =  0 . 0  ; 

W~b  nrml  [0]  :=  W  b  nrml  [1]  ; 
vTahray  [o]  :=  w  Array  [i]  ; 

THERM_OON_ARRAY  To]  :=  THERM_CCN_ARRAY  [1]  ; 

(*  checking  for  w_array  values  of  zero  so  as  to  set 

u_a,  u  v,  h_a_v,  h_v_a  to  zero.  ALL  w_array  values 

must”be  zero  for  u_a,  u_v,  h_a_v,  h_v_a  to  be  changed  *) 

W_zero_flag  :=  true; 

FOR  I  ;=  1  TO  N  DO 

IF  W_ARRAY  [I]  >  0  THEN 

W_zero_flag  ;=  false; 

K_AVG  :=  0.0; 

FOR  I  :=  0  TO  N-l  DO 

K_AVG  :=  K  AVG  +  THERM  CCN  ARRAY  [I]  * 

(sqr(RAdial_dist  Array  [i+i])  -  sqr (radial  dist_array  [i])),- 


(*  N.B.  . . - . - - - - - - - 

re-assign  w_array  to  W_TISSUE 

REASON  :  necessity  to  add  n-l  tissue  layers  before  interface 
adjustments  have  been  made.  In  preparation  for  calculation 
of  capillary  data.  Values  are  found  at  exact  tissues  locations, 
not  one  value  greater  than  actual  tissue  location. 

- - — - - - — . - . --*) 

for  i  :=  1  to  N-l  do 

WJTISSUE  [I]  :=  w_array  [i+1]  ; 
w_tissue[0]  :=  w_array  [1] ; 

(*  apply  equalizing  formula  on  interface  locations 
p  =  {GAMBIA  [i]  *  p-  +  p)  /  (1  +  GAMBIA  [i] ) 
also  determine  B_i_l  and  ALPHA_PROP  arrays  *) 

FOR  I  :=  1  TO  N  DO 

IF  I  IN  I  INTERFACES  THEN 
BEGIN 

GAMM  :=  GANNA  (I,H  r_ARRAY)  ; 

1HERM_CCN  ARRAY  [IT  :=  (GANN  *  THERM  CCN  ARRAY  [1-1]  + 

THERM  CCN  ARRAY  Tl+lJ)  /  (1  +  GANN)  ; 
HEAT_CAP_ARRAY  [I]  :=  (2ANN  *  HEAT  CAP  ARRAY  [1-1]  + 

HEAT  CAP  ARRAY  [T+1]T  /  (1  +  GANN)  ; 
DENSITY  ARRAY  [I]  :=  (SAW!-*  DENSITY  ARRAY  [1-1]  + 
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DENSITY  ARRAY  [1+1] )  /  (1  +  GAT-M)  ; 

W  ARRAY  [I]  (GANM  *  W  AHRAY[I-l]  +  W  ARRAY [1+1])  /  (1+GAhW)  ; 
VTb  nrml  [I]  :=  W_ARRAY[T]  *  W  B  NORMALIZING  VAL; 

Q_XRRAY_PHYS  [I]  :=  (GANN  *  Q  ARR5Y  PHYS  [I-lJ  + 

Q_ARRAlTPHYS  Tl+1] )  /  (1  +  GARM)  ; 

END; 


FOR  I  :■  1  TO  N  DO 
BEGIN 

B_i_l  [I]  :=  (H_l  *  LEN)  /  1HERM_OON_ARRAY  [I]; 

ALPHA_PROP  [I]  :=  THERM  CCN_ARRAY  [I]  /  (DENSITY  ARRAY  [I]  * 

Heat  cap  array  [i] )  ,• 

Q  ARRAY_nrml  [I]  :•=  Q  ARRAY“PHY5  [I]  *  SQR(PHYS  RAD)  / 
(NORMALIZING  TEMP  *  BLOOD  THERM  CON)  ; 

(*  :=  0.0;  *)  -  - 

END; 


(*  establish  U_a  and  U_v  +) 

IF  W  zero_flag  THEN 
BEGIN 

U  A  :  =  0.0; 
tTV  :=  0.0; 

H~A  V  :=  0.0; 

H_\TA  :=  0.0; 

END 

ELSE 

BEGIN 

u_a_phys  :=  2  *  k_avg  /  sqr  (phys_rad)  /  In  (2  *  sqrt  (r_v/r_a) )  ; 
U_v_j)hys  :=  2  *  k_avg  /  sqr  (phys_rad)  /  In  (2  *  sqrt  (r_a/r_v) ) ; 

U  a  :=  (U_aj?hys  *  SQR  (PHYS  RAD)  /  ORG  VALS  [BBLOOD]  .THERM  OCN)  ; 
l Tv  :=  (U_v_phys  *  SQR  (PHYSJRAD)  /  ORGJVALS  [BBLOOD]  .THERM_CCN)  ; 

U  A  DIV  :=  u  a  /  (n+m)  ; 

CnTDIV  :=  tTV  /  (n+m)  ; 

(*  establish  H_a_v  *) 

H_a_v_PHYS  :=  pi*  phys  h  z  *  2  *  K  AVG  / 

HYPERBOLJOOS  T7  -  (SQR (R  V  -  R  a)  /  (2  *  R  a  *  R  v) ) ) ; 
H  A_V  :=  H  A  V  PHYS  /  (PI  *  ORG  VALS [BBLOOD]  .IHERM  CON  *  ~ 

RifS  H  Z)  *  A  a* 

H  V  A  :=  H  A  V PHYS  /  TPI  *  ORG  VALS  [BBLOOD]  .IHERM  CON  * 

FRY3_H_Z) *  A_v; 

END; 

{  reset  u_a,  u_v,  h_a_v,  h_v_a, } 

(* 

h_a_v  :=  0.0; 
h  v  a  ;=  0.0; 

*7 


(*  U  V  DIVIDED  BY  5  TO  ACCOUNT  FOR  COUNTERCURRENT  HEAT  EXCHANGE  BETWEEN 
THE  BLOOD  VESSELS  *) 

IF  (H  V  A  <>  0)  AND  (H_A  V  <>  0)  THEN 
BEGIN- 

U  V  :=  U  V/5.0; 

UTV_DIV  :=  U_V_DIV  /  5.0; 
end; 


(*  establish  m_a  and  m_v  arrays  *) 
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(*  establish  DIMENSIONAL  capillary  loss  *) 

CAPIL_SUM  :=  0.0; 

PC®  I  :=  0  TO  n-1  DO 

CAPIL  SUM  :=  CAPIL  SUM  +  W  TISSUE  [I]  * 

(SQR  (PHYS"R_i  [I+lT)  -  SQR(PHYS_P_i[I] ) )  ; 

CAPILLARY  LOSS  PI  *  PHYS  Hz*  CAPIL  SUM; 

START  M  a  .-^CAPILLARY  LOSS  T  "AXIAL  DlVS; 

M_a[lT  :=  START_M_a  ; 

FOR  J  :=  2  TO  M  DO 

M_a  [J]  :«  M_a  [J-l]  -  CAPILLARY_LOSS ; 

(*  controlling  the  negative  zero  *) 

IF  M  a  [M]  <  CAPILLARY  LOSS  THEN 
M>  [M]  :=  0.0; 

M  V  [M]  :=  0; 

FOR  J  :=  M-l  DCWNTO  1  DO 

M_v  [J]  :=  M_v  [J+l]  +  CAPILLARY_L0SS; 

(*  establish  B_a  and  B_v  arrays  *) 

FOR  J  :=  1  TO  M  DO 
BEGIN 

B  a  [J]  :=  (M  a  [J]  *  SQR  (PHYS  RAD))  / 

(Org  vals  [bbloqd]  .Density  *  pi  *  sqr  (r  a)  * 

PHYS  H  Z  *  BL  ALPHA) ; 

B  V  [J]  :=  (M  v  [Jl  *  SQR  (5 HYS  RAD))  / 

(OKG  VALS  [BBLOQD] .DENSITY  *  PI  *  SQR  (R  v)  * 
P0YS_H_Z  *  BL_ALPHA) ; 

END; 

(*  establish  B  i 

N.B.  -  B_T_C  in  data  dump  routine  is  B_i  in  the  program  code  *) 
B_i  :=  (h_c  *  PHYS_RAD)  /  ORGJVALS  [SSKIN]  .THERM_OON; 

(*b_i  :=  0.0;*) 


(*  calculate  dimensionless  DELTA  *) 

SUM  :=  0.0; 

FOR  I  :=  0  TO  N-1  DO 

SUM  :=  SUM  +  W  tissue  [I]  *  (SQR (RADIAL  DIST_ARRAY[I+1] )  - 
SQR (RADIAL_DIST_ARRAY [I] )T; 

total_a  :=  -2  *  b_a  [1]  +  A_a  *  SUM  -  (A_A  *  U_a  +  H_a_v) ; 

total_v  :=  -2  *  b  v  [1]  -  2  *  A_v  *  SUM  -  (A_V  *  U_v  +  H_v_a)  ; 
total  :=  abs  (totaI_v) ; 

if  abs  (total  a)  >  abs  (total_v)  then 
total  :=  abs  (total_a) ; 

DELTA  :=  0.4  /  TOTAL  ; 


fdeltaset} 

{max  delta  not  to  exceed  0.004} 

MAX  DELTA  :=  DELTA; 

DELTA_SEC  :=  GC_DELTA; 

WRITElN; 

WRITELN; 

WRITELN  ( ' ORIGINAL  DELTA  IN  SECONDS 1 , DELTA  SEC:  10: 6)  ; 

WRITELN; 

WRITELN  ( 'PROGRAM  CALCULATED  MAX  DELTA  VALUE  =  ' ,  MAX  DELTA  :  20:15)  ; 
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WRITELN; 

DELTA  :>  DELTA  SEC  *  BL  ALPHA  /  SQRfPHYS  RAD)  ; 
writeln  (’USER  CALCULATED  delta  =  ' ,DELlS:20;15)  ; 
WRITELN; 

PAUSE  CRT; 


(*  loading  numeric  constants  into  variables  of  type  EXTENDED  *) 

one  :=  c_one; 
two  :*  c_two ; 
three  :  =  c_three ; 
eight  :*  c_eight; 

(*  normalize  w  tissue  for  use  in  t_a_t_v  calculations 

which  require  an  absence  of  gairma  values  used  in  interface 
locations  *) 

FOR  I  :=  0  TO  N-l  DO 

W_TISSUE [I]  :=  W_TISSUE  [I]  *  W_B JJORMALI  ZIN3_VAL ; 

END;  (*initialize  *) 


BEGIN  (*  main  program  *) 
clrscr; 

WRITE  ( 'HOW  MANY  TIMESTEPS  ?  ' )  ; 

READLN  (LCV_R) ; 

WRITELN; 

WRITELN ; 

WRITE  ( '  intervals  at  which  to  print  tenperature  data  ?  ' )  ; 

READLN  (INTERVAL)  ; 

WRITELN; 

ANATOMICAL_INPUT  (DIAM,  LEN)  ; 

NUM_DIVISICNS  (OOREJ3IVS,  MUSCLEJDIVS,  FAT_DIVS,  SKIN_DIVS,  N,  M)  ; 

WIDTH  CALCULATIONS  (H  r  CORE,H_r  MUS , H_r_FAT , H  rjSKIN,  CORE  MUS_GANMA, 

~  MUS  EAT  GAMMA,  EAT  SKIN  GAiWA,  DIAM,  LESJ, 

CORE_DI VS,  MUSCLE J5IVS , FAT_DIVS ,  SKIN_DIVS) ; 

EST  INTERFACES  (I  INTERFACES,  I  CORE  MUS,  I  MUS  FAT,  I  FAT  SKIN, 

SKINJSUREACE,  CDREJDTVS,  MU5CLET>IV5,  EATJDIVS,  SKIN_DIVS)  ; 

CREATE  ARRAYS  (H  r  ARRAY,  RADIAL  DIST  ARRAY,  I  CORE  MUS,  I_MUS  FAT, 

I  EAT  SKIN,  H_r_CORE,  H_r_MUS,  H_r_FAT,  H_r_SKUJ) ; 


DESTINATION; 

INITIALIZE; 

TEMP  INITIALIZE  ; 


DATA  DUMP; 


TISSUE; 

CLOSE  (F) ; 

CLOSE  (F2) ; 

CLOSE  (COT_DATA)  ; 

CLRSCR; 

WRITE  ('PROGRAM  COMPLETED - '); 

PADSE_CRT; 

END.  (*  main  program  *) 
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2.  VAR  SP.PAS 


{GLOBAL  DECLARATION  UNIT 


INCLUDES  CONST,  TYPE  AND  VAR  DECLARATIONS  } 


UNIT  varjsp; 

INTERFACE 
Uses  CRT; 

CONST 

{///////////////////  begin  control  point  section  \\\\\\\\\\\\\\\\\\\\\\\\\\\ } 

(*  GC_DELTA  =0.0  will  allow  program  to  determine  delta  value.* 

GC  DELTA  set  to  any  other  value  will  over-ride  program 
calculation  of  delta  value. 

N.B.  -  GC  DELTA  MUST  BE  GIVEN  A  VALUE  IF  ALL  W  VALUES  HAVE  BEEN 

Set  to  zero.  *) 

(*  DELTA  IS  A  FULL  TIME  STEP  -  it  will  be  adjusted  later  in 
the  tissue  equations  to  be  a  half  time  step.  In  the 
blood  equations  DELTA  will  be  a  half  time  step  for  the  first 
iteration  of  the  loop  and  a  full  time  step  on  all  subsequent 
iterations . 


N.B.  --  GC_DELTA  REPRESENTS  DELTA  IN  SECONDS  *) 

GC_DELTA  =1.0;  (*  PROGRAM  MINIMUM  =  0.0000015;  *) 

(*  length  and  diameter  must  be  in  CENTIMETERS  *) 

GC  LEN  =  8.0; 

GCJDIA  =  1.5; 

(*  control  point  for  number  of  radial  divisions  *) 

GC_CORE_D  TVS  =  3; 

GC_MUS  DIVS  =  3; 

GC_FAT  DIVS  =  2; 

GC_SKIN_DIVS  =  4; 

(*  control  point  for  the  number  of  axial  divisions  *) 

GC_AXIAL_DIVS  =  10; 

(*  control  point  for  temperature  values  in  degrees  C  *) 

TA  =  35;  (*  initial  assignment  for  j=l  to  j_max-l  values  *) 

(*  ta  at  j=l  held  constant  during  program  run  at  initial  value*) 

TV  =  35;  (*  used  for  initial  assignment  at  j  max-2  to  j  =  1  values  *) 

(*  seed  value  at  j_max-l  is  calculated  by  program  *) 

MESH_TEMP  =  35; (*  used  to  set  initial  tissue  temperature  throughout  mesh*) 

NORMALI Z ING_TEMP  =  37; 

hi  =  8.09;  (*  units  =  Watt  /  meter  squared  /  deg.  C  *) 
h  C  =  5.17;  (*7.02;*) 
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T  0  = 


0 .0/normalizing_tenp  ;  (*  environmental  temperature  ,  deg.  C  *) 


(*  basal  metabolic  rate  -  CONTROL  POINT  FOR  Q  units  =  Watt  /  cu.m  *) 

core_metab  =  170.5; 
musjmetab  =  631.9; 
f at_metab  =  5.0; 
skin  rnetab  =  247.4; 


(*  basal  blood  flow  rate  CONTROL  POINT  for  W  units  =  kg  /  sec  /  cu.m 
N.B.  -  if  ALL  of  these  values  are  set  to  zero  then  GC_DELTA  MUST  BE 
ASSIGNED  A  VALUE  at  the  top  of  this  section.  *) 


core  blood_meters  =  0.173; 
mus  El ood  meters  =  0.641; 
fat^bloodT meters  =  0.0; 

(*  to  account  for  a  zero  blood  flow  at  the  i=n  value  *) 

skin_blood_meters  =  0.251  *  gc_skin_divs/  (gc_skin_divs-l)  ; 

(*  organ  radius  to  finger  ratio  (R_i  /  R)  *) 

core  ratio  =  0.7057; 
muscTe_ratio  =  0.7954; 
fat_ratio  =  0.8099; 
skin_ratio  =  1.0; 

(*  thermal  conductivity  units  =  Watt  /  m  /  deg.C  *) 

core  therm_con  =  1.064; 
mus_therm_con  =  0.418; 
f at_therm_con  =  0 . 204 ; 
skin  therm_con  =  0.293; 
blood_therm_con  =  0.45; 

(*  heat  capacity  units  =  J  /  kg  /  deg.  C  *) 


core  heat_cap  =  2102.0; 
mus  Eeat_cap  =  3136.0; 
fat~heat_cap  =  2520 . 0 ; 
skin  heat_cap  =  3780.0; 
bloo3_heat_cap  =  3899.0; 

(*  density  units  =  kg  /  cu.  m  *} 

core  density  =  1401.0; 
musjdensity  =  1057 . 0 ; 

f at_density  =  900 . 0 ; 
skin  density  =  1057.0; 
bloo3_density  =  1060.0; 

(*  Radii  of  artery  and  vein  -  units  =  meters  *) 


R_A  =  1000E-6; 
R  V  =  1500E-6; 


/////////////////////  end  control  point  section  \\\\\\\\\\\\\\\\\\\\\\\\\\\ 
//////////////  NOTHING  should  be  altered  below  this  point  \\\\\\\\\\\\\\\\ 


J_MAX  =  GC_AXIAL  DIVS  +  1; 

I  MAX  =  GC  COREJ5IVS  +  GC_MUS_DIVS  +  GC_FAT_DIVS  +  GC_SKIN_DIVS  +  1; 
MSX_PTS  =  Tjmax  *  j_max  -  i_max; 

(*  values  to  be  converted  to  extended  precisions  variables  later  *) 
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C_two  =  2.0; 
c_eight  =  8.0; 
c_three  =  3.0; 
c_one  =  1.0; 

TYPE 

PTR  TYPE  =  "TCM  NODE; 

TOM- NODE  =  RECORD 
TNFO  :  EXTENDED; 

NEXT,  BACK  :  PTR_TYPE; 

END; 

GRID  =  ARRAY  [0.  .I_MAX,  1..J  MAX]  OF  extended; 

BOUNDS  =  SET  OF  1. .1  MAX; 

extended  ARRAY  =  ARRSY  [1.  .MAX_PTS]  OF  extended; 

J  ARRAY  =  ARRAY  [1..J  MAX]  OF  extended; 

I_ARRAY  =  ARRAY  [0..I  MAX]  OF  extended; 

A_TYPE  =  (A_FIRST,  A  PEG,  A  END,  A  EXTERNAL,  A  EXTERNALJ2ND)  ; 

R_TYPE  =  (R_CENTR  El®,  R_CEf?TER_RE5 ,  R_EXTERNAL_END , 

R_EXTERKAIjREG(  R_INTERFACE  REG,  R  INTERFACE_END ,  R  REG_END) ; 
S  TYPE  =  (S  REG_CENTER_INT ,  S_EXTERNAD,  S_E ND,  S_END_E$TERNALT  ; 

CO  ORD  =  RECORD 

I_LOC,  J_LOC  :  INTEGER; 
tom_tenp  :  extended; 

END; 

L  TYPE  =  ARRAY  [l..MAX_PTS]  OF  CO_ORD; 

COEF_FIELDS  =  RECORD 

MINUS,  STD,  PLUS  :  extended; 

END; 

COEF  ARRAY  =  ARRAY  [1 .  .  I_MAX,  1 .  .  J_MAX]  OF  COEFJFIELDS; 

VALUES  =  RECORD 

THERM  CON,  HEAT  CAP,  DENSITY, 

METAB,  METERS_BSS_BLOOD  : extended; 

END; 

ORGAN  =  (CCORE,  MMUSCLE,  FFAT,  SSKIN,  BBLOOD)  ; 

VAL_ARRAY  =  ARRAY  [ORGAN]  OF  VALUES; 

STR  TYPE  =  STRING  [20]  ; 


VAR 

U_PTR,  UTRV,  P,  PI  :  PTR_TYPE; 
av_step,  half  av  step,t_step,  LCV_R,  INTERVAL  :real; 

N  ,  M,  AXIAL  DIVE,  ROW  TIME,  I,  J, 

OORE_DIVS ,  MOSCLEJDIVS ,  FAT  DIVS, 

SKIN_DIVS,  SKIN_SURFACE,  D  K,  test, lev, 

count,  I_CORE_MUS ,  I_MUS_F5T,  I_FAT_SKIN,  COEF_SUB  :  INTEGER; 

DIAM,  LEN,  H_Z ,  PHYS  H_Z,  RAD,  PHYS_RAD, 

PHYS  LEN,  U  a,  U_v,  H  a_v,  H_V_A,  H_A_V_PHYS, 

DELTA,  MAX_DELTA,  DELl!A_sec,  A_a,  A_v, 

BL  ALPHA,  W  b  normalizing  val, 

B  T,  CAPILLARY  LOSS,  U_a_phys,  U_v_phys, 

START  M  A,  CAPTL  SUM, 

H  r  CORE,  H_r  MUS,  H_r_FAT,  H_r_SKIN, 

CORE  MUS  GAMMA,  MUS  FAT  GAMMA,  FAT  SKIN  GAMMA, 

core_widTh,  mlt_wieTh,  Eat  width,  EkinJJidth, 

CORE  MUS  BOUND,  MUS  FAT  BOOND, 
fat_Ekin~bound,  SKlfl  SURFACE_BOUND,  k_avg, 

one,  two,  three,  eigHt,  zero,u_a_div,  u_v_div  :  extended; 

I_INTERFACES,  CORE_SET,  MUS_SET,  FAT_SET,  SKIN_SET  :  BOUNDS; 

RADIAL  DIST  ARRAY,  H  r  ARRAY,  W  ARRAY,  W  TISSUE, 

PHYS_R  i,  PHYS  H_r  ARRAY,  W  b  nrml,  Q  AREAY_nrml, 

THERM  0OJ  ARRAY,  HEAT  CAP_ARrSY,  Q_ARRAY_PHYS , 

DENSITY_AERAY,  ALPHA_PROP,  B_i_l  :  I_ARRAY; 
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A  NODE 
RHODE 
S~NQDE 


A_TYPE 
R_TYPE 
S  TYPE 


B  a,  B  v,  M  a,  M  v, 

Try.  T"A,  TTAjoar,  T  V  bar, 

t_a_no3e ,  t_v_node ,  ol3_t_a_node ,  old_t_v_node  :  J_ARRAY; 

TEMP  GRID,  HALF_GRID  :  GRID; 


LOCATE  :  L_TYPE; 

f , f 2 , OOT_DATA  :  TEXT; 


ORG  VALS 


VAL  ARRAY; 


RESP2 , GRAPH_RESP  :  CHAR; 

tom_flag,w_zero_flag  :  boolean; 


FUNCTION  GANMA  (I  VAL:  INTEGER;  VAR  H_r  ARRAY  :  I  ARRAY)  :  extended; 
PROCEDURE  PAUSE  CRT; 

PROCEDURE  DERIVE  CENTER  TEMP  {VAR  DEGREES  .-GRID)  ; 

PROCEDURE  DUMP  KELP  ARREYS ; 

PROCEDURE  DUMP  T  ARRAY  (STEP  VAL  :  REAL;  DEGREES  :  GRID)  ; 

PROCEDURE  CALCTTS_TV  half  (STEP  :  EXTENDED)  ; 


IMPLEMENTATICN 


FUNCTICN  GAMMA  (I_VAL:  INTEGER;  VAR  H_r_ARRAY  :  I_ARRAY)  :  extended; 

{ . returns  the  GANMA  value  for  any  I  location  in  the  mesh . 

BEGIN 


IF  I  VAL  >=  N  THEN 
GEMMA  :=  0.0 
ELSE 

IF  I  VAL  =  1  THEN 
GEMMA  :=  1.0 
ELSE 

GAMMA  :=  H_r_ARRAY [ I_VAL]  /  H_r_ARRAY [I_VAL  +  1] ; 

END; 


PROCEDURE  PAUSE_CRT; 

{ . halts  screen  output  scrolling 

VAR 

AGAIN  :  CHAR; 

BEGIN 

WRITE  ( '  PRESS  ANY  KEY  TO  CONTINUE  :  ' )  ; 

AGAIN  :=  READKEY; 

WRITELN; 

CLRSCR; 

END; 


PROCEDURE  DERTVE_CENTER_TEMP  (VAR  DEGREES  :GRID) ; 
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{  establishes  center  temperatures  at  t  [0 , j ] 

*****  N.B.  this  procedure  REQUIRES  A  MINIMUM  OF  2  RADIAL  DIVISIONS  } 

VAR 

J,  I  :  INTEGER; 

BEGIN 

FOR  J  :=  1  TO  M  DO 

DEGREES  [0,Jj  :=  DEGREES  [2,J]  +  1.5  * 

(DEGREES  [1,J]  -  DEGREES  [2,J]); 

END; 


PROCEDURE  DUMP_HALF_ARRAYS; 

{ . prints  t_a,  t_v  arrays . } 

VAR 

J_val  :  INTEGER ; 

BEGIN 

CLRSCR; 

WRITELN  (F, 'TIMESTEP  =  HALF_AV_STEP  :  5:2,  ' 

★  *  )  • 

WRITE  (F, 'T_A' ) ; 

FOR  J  VAL  :=  1  TO  M-l  DO 

WRITE  (F,  J_VAL:3,  ')  1  ,T_A  NODE  [ J  VAL]  *  NORMALIZING  TEMP:7:3)  ; 
WRITELN (F) ; 

WRITELN (F) ; 

WRITE  (F, 'T_V' ) ; 

for  J  VAL  1  to  m-l  do 

write  ( f , j_val : 3 , ' ) ' , t_V_node [ j_val]  *  NORMALIZING_TEMP :  7 : 3 )  ; 
writeln(f ) ; 

END; 


PROCEDURE  DUMP_T_ARRAY  (STEP_VAL  :  REAL;  DEGREES  :  GRID) ; 

{ . . .prints  two  dimensional  temperature  matrix  or  half_step  matrix. . } 

VAR 

I,J  :  INTEGER; 

BEGIN 

WRITELN  (F)  ; 

WRITELN  (F, 'TEMPERATURE  ARRAY  AT  TIMESTEP  '  ,  STEP  VAL:7:2,  ’  (end/ext)  ' :42) ; 
FOR  I  :=  N  DOWNTO  1  DO 
BEGIN 

(*  WRITE  (F,  ' 1= ' , I) ;  *) 

FOR  J  :=  1  TO  M  DO 

WRITE  (F,  DEGREES  [I,J]  *  normalizing  temp:7:3); 

WRITELN  (F)  ; 

END; 

WRITELN  (F) ; 
pause_crt ; 

END; 


PROCEDURE  CALC_TA_TV_half  (STEP  :  EXTENDED) ; 

{....calculates  ta_bar[j]  and  tv  bar[j]  in  preparation  for 

s_vector  calculation  . } 
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:  extended; 


VAR 

SUMl,  SUM2,  SUM3,W_VAL.  HHaa  w, 

T  v  plus,  B  v  plus,  a,  b,  TEMP,  DELTA_VAR 
iT^Iyal  ~  :  INTEGER; 

first, second, third, answer  :  extended; 

BEGIN 

IF  STEP  =0.0  THEN 

DELTA_VAR  :=  DELTA  /  2 
ELSE 

DELTA_VAR  :=  DELTA; 

FOR  J  VAL  :=  2  TO  M-l  DO 
BEGIN 

(*  CALCULATE  T_A_BAR  *) 

SUMl  :=  0.0; 

SUM2  :=  0.0; 

SUM3  :=  0.0; 

FOR  I  :=  0  TO  N-l  DO 
BEGIN 

W_VAL  :=  W_TISSUE  [I]; 

SUMl  :=  SUMl  +  W  val  *  {SQR  (RADIAL  DIST  ARRAY [1+1] )  - 

SQR  ( RAD IAL_DI ST  ARRAY  [I] ) ) ; 

SUM2  :=  SUM2  +  (TEMP  GRID[I,J  VAL-lJ  +  TEMP  GRID  [I,J  VAL])  * 
(SQR  (RADIAL  Dl5T_ARRAY [1+1] )  - 
SQR  (RADIALJDIST_ARRAY [I] ) ) ; 


SUM3  :=  SUM3  +  W  VAL  *  (TEMP  GRID  [I,  J  VAL-1]  + 

TEMP_GRlD  [I,J  VAL]T  *  (SQR  (R5DIAL_DIST_ARRAY  [1+1] ) 
-  SQR  (RADIAL  DIST_ARRAY  [I] ) ) ; 


END; 


(*  FINAL  FORMULAS  *) 

(*  T_A  is  the  t_a_BAR  value  at  the  previous  full  time  step 
the  same  is  true  with  regard  to  T_V  and  t_v_BAR  *) 

T_a_bar  [J_VAL]  :=  t_a  [j_val]  +  DELTA_V7^R  *  ((  -two  *  B_A  [J_VAL-1] 

A_A  *  SUMl  -  A_A  *  U_A  -  H_A  V)  *  T_A  [J_VAL]  + 

(two  *  B  A  [J_VAL-1]  -  A_A  *  SUMl)  *  T_A  node  [JVAL-1] 
+0.5*  5_A  *  U_A  *  SUM2  +  H_A_V  *  T_V  [J_VAL] ) ; 

END;  (*  J-LOOP  *) 


(*  CALCULATE  T_V_BAR  *) 

FOR  J  VAL  ;=  2  TO  M-l  DO 
BEGIN 

SUMl  :=  0.0; 

SUM2  :=  0.0; 

SUM3  ;=  0.0; 

FOR  I  :=  0  TO  N-l  DO 
BEGIN 

W_VAL  :=  W_TISSUE  [I]  ; 

SUMl  :=  SUMl  +  W_val  *  (SQR  ( RADIAL_DI ST_ARRAY [ I + 1 ] )  - 
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SQR  {RADIAL  DIST  ARRAY  [I] ) ) ; 

SUM2  :=  SUM2  +  (TEMP  GRID[I,J  VAL+lJ  +  TEMP  GRID  [I,J  VAL])  * 

{SQR  (Radial  diSt  array [i+i])  - 

SQR  ( RADIAL J5l  ST  ARRAY  [I]  ) )  ; 


SUM3  :=  SUM3  +  W  VAL  *  (TEMP  GRID  [I,  J  VAL+1]  + 
TEMP  GRID  [I, J_VAL] )  * 

(S§R  (RADIAL  DIST  ARRAY  [1+1] ) 

-  SQR  (RADIAL  DIST  ARRAY  [I] ) )  ; 

END;  (*  I-LOOP  *) 


T_v_bar [ J_VAL]  :=  t  v  [j  val]  +  DELTA  VAR  *  ((two  *  B  V  (J  VAL-l] 

-  two  *  A  \f~*  SUMl  -  A _V  *  U_\7  -  H  V  A)  *  T  V  [J  VAL]  + 

(-two  *  E  V  [J_VAL-1]  +  A _V  *  SUMl)-*  T  V  RODE  TJ  VAL-l]  + 
0;5*AV*UV*  STJM2  +  H  V  A  *  T  A  [3  ^AL]  +  075  *  A  V  * 
SUM3)  ;  “ 


END;  (*J  VAL*) 


(*  calculate  t_a_node  and  t_v_node  *) 
for  j_val  : =  2  to  m-1  do 

t_a_node  [j_val]  :=  2  *  t_a_bar  [j_val]  -  t_a_node  [j_val-l] ; 


SUMl  :=  0.0; 

FOR  J_VAL  :  =  M-1  TO  M  DO 
FOR  I  :=  1  TO  N  DO 

SUMl  :=  SUMl  +  HALF  GRID  [I, J_VAL] ; 
T  V  NODE  [M-1]  :=  SUMl  /  (2*  N)  ; 


for  j_val  :=  m-1  downto  2  do 

t_v_node  [j_val-l]  :=  2  *  t_v_bar [ j_val]  -  t_v_node  [j_val] ; 

END;  (*  calc_ta_tv_half  *) 

END . 
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3.  TISUE  SP.PAS 


Uhit  tisue_sp; 

interface 

PROCEDURE  TISSUE; 

implementation 

uses  var_sp,  phys_sp,  crt,  dos,  printer; 


PROCEDURE  TISSUE; 

VSR 

ORIG  A  r  COEF,  ORIG  A  z  COEF  :  COEF  ARRAY; 


PROCEDURE  AXIAL  DIR  (  I  VAL,  J  VAL  :  INTEGER; 

VAR  J_MENUS_PORM ,  J_FORM,  J_PLUS_PORM  :  extended)  ; 

{ . imports  the  I  and  J  coordinates,  identifies  the  type 

of  node  using  enumerated  data  types  and  calculates  three 
coefficients  at  the  given  I,J  location  to  fill  the 
ORIG  A  z  COEF  array  used  for  the  traversal  in  the 
RADIAL  direction . } 


VAR 

HEAT_DEN_FRACT,  AjSQRD,  FRACT,  W_U_A,  tenp_w_u_a  :  extended; 
BEGIN 

(*  establish  node  type  *) 


IF  J  VAL  =  1  THEN 

5  NODE  :=  A  FIRST ; 

IF  J  VAL  =  M  THEN 
IF  I  VAL  <>  N  THEN 

a  Rode  :=  a_end 

A  NODE  :=  A  EXTERNAL  END; 
IF  (J- VAL  <>  M)  AND  (J  VEL  <> 
IF~I_VAL  o  N  THEN  ~ 
A_NCDE  :=  A  REG 
ELSE 

A  NODE  :=  A  EXTERNAL; 


{pt. 

0} 

{pt. 

4,6,7} 

{pt. 

5} 

THEN 

{pt. 

1-3} 

{pt. 

2} 

(*  preliminary  calculations  *} 


A  SQRD  :=  (BL  ALPHA  /  ALPHA  PROP  [I  VAL] )  *  SQR (LEN/PHYS  RAD)  ; 

HEAT  DEN  FRACT  :=  (ORG  VALSlBBIOOD] 7DENSITY  *  ORG  VALS  [BELOOJ]  .HEAT  CAP)  / 

TDENSITY  ARRAY  [I  VAL]  *  HEAT  CAP  ARRAY  [I  VALT)  ; 
FRACT  :=  one/ (A  SQRD  *  SQR(H  zT)  ; 


W_U_A  :=  W_b_nrml  [I_VAL]  +  U_A_div  +  U_v_div; 

(+  setting  w,  u_a,  and  u_v  at  external  skin  level  to  0.0  *) 
tenp_w_u_a  :=  0.0; 

(*  FINAL  CALCULATICNS  *) 

CASE  A  NODE  OF 
A  FIRST  :  ; 


A  REG  :  BEGIN 

J  MINUS_FORM  :=  FRACT; 

CT FORM  :=  -one*  (two  *  FRACT)  -  (HEAT  DEN  FRACT  *  (W  U  A  /  two) )  ; 
OTPLUS_FORM  :=  FRACT;  “  ” 

END; 
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A  EXTERNAL  :  BEGIN 

J  MINUS  FORM  :=  FRACT; 

J"~ FORM  :=  -two  *  FRACT  -  (HEAT  DEN  FRACT  * 

(three  *  tenp  W  U  A)  /  eight]"; 

J_PLUS_FORM  :=  FKAW7 

END; 

A  END  :  BEGIN 

J  MINUS_FORM  :*  two  *  FRACT; 

FORM  :=  -two  *  FRACT  *  (one  +  H  z  *  B  i  1  [I  VAL) ) 

-  ( (HEAT  DEN  FRACT  T  W  b  nrml  [I~VAL] )  /  two  )  ; 
J_PLUS_FORM  :=  0.0; 

END; 

A  EXTERNAL  END  :  BEGIN 

J  MINUS  FORM  :=  two  *  FRACT; 

CT FORM  :=  -two  *  FRACT  *  (one  +  H  Z  *  B  i  1  [I  VAL] )  - 
~  (three  *  HEAT  DEN  FRACT  *  W  b  nrml  [T  7AL]  “/  eight) ; 
J  PLUS_P0RM  :=  0.0; 

EKD; 

END;  (*CASE*) 

END;  (*  axial__dir  *) 

PROCEDURE  RAD  DIR  (I  VAL,  J  VAL;  INTEGER; 

VAR  T_MINUS_P0RM,  I_F0RM,  I_PLUS_FORM  :  extended)  ; 

{ . inports  the  I  and  J  coordinates,  identifies  the 

type  of  node  using  enumerated  data  types  and  calculates 
three  coefficients  at  the  given  I,J  location  to  fill  the 
ORIG  A  r  COEF  array  used  for  the  traversal  in  the 
SECX^~TRAVERSAL . } 

VAR 

ALPHA  FORM,  HEAT  DEN  FRACT, 

FRACTT  FORM,  GAT*!,  IQ?_I,  R_I ,  tenp_f orm  ;  extended; 

BEGIN 

(*  establish  node  type  *) 

IF  (J  VAL  =  M)  THEN 
BEGIN 

R  NODE  :=  R  REG  ENE,  {pt.  4} 

IP  (I  VAL  HI  ^INTERFACES)  THEN 

r  Node  :=  r  interface  end,-  {pt.  7} 

IF  (I  VAL  =  1)~THEN 

R  RODE  :=  R  CENIR  END;  {pt.  6} 

IF  I~VAL  =  N  THEN 

RHODE  :=  R  EXTERNAL  END;  {pt.  5} 

END 

BEGIN 

IF  I  VAL  =  1  THEN 

R  Rode  :=  R  CENTER  REG;  {pt.  1} 

IF  T  VAL  =  N  THEN  ” 

R  RODE  :=  R  EXTERNAL  REG;  {pt.  2} 

*  Tl  VAL  <>  T)  AND  (I""'VAL  <>  N)  THEN 
R_NGDE  ;=  R_INTERFACS_REG;  {pt.  3} 

END; 

(*  preliminary  calculations  *) 

ALPHA  FORM  ALPHA  PROP  [I  VAL]  /  BL  ALPHA  ; 

HEATJDEN  FRACT  :=  (GkG  VALSlBBLOOD]  .DENSITY  * 

ORG  VALS [BBIOOD] .HEAT  CAP)  / 

(DENS ITYjERRAY [ I_VAL]  *  HE5T_CAP_ARRAY[I_VAL] ) ; 

(*  H_R  ARRAY  AT  THE  1  LOCATICN  MUST  USE  THE  VALUE  AT  H_R_ARRAY  [2]  *) 

IF  I  VAL  =  1  THEN 
BEGIN 

FRACT  :=  two  /  SQR(H  R  ARRAY  [2] )  ; 
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GAM-i  GAM1A  ( 2 ,  H  R  ARRAY)  ; 

H_R_I  :=  H_R_ARRAY  j2j ; 

END 

KI 

BEGIN 

BRACT  :*  two  /  SQR(H  R  ARRAY  [I  VAL] ) ; 
GAM»I  :=  GAM1A  (I  VAL,  H  R  ARRAY)  ; 
H_R_I  :=  H_R_ARR5Y  [I_V5LT; 

END;  ~ 


FORM  :=  W_b_nrml [I_VAL]  +  U_a_div  +  U_v_div; 

(*  setting  w,  u  a,  and  u_v  at  external  skin  level  to  0.0  *) 
terrp_form  :*  O.TS; 


R_I  :*  RADIAL_DIST_ARRAY  [I_VAL]  ; 

(*  final  calculations  *) 

CASE  R  NODE  OF 

R  CENTER  REG  :  BEGIN 

“  I_MINUS_FORM  :=  0.0; 

INFORM  :=  -one  *  ALPHA  FORM  *  FRACT  -  (HEAT_DEN  FRACT 
*  FORM  /  two)  ; 

I  PLUS_P0RM  :=  ALPHA_PORM  *  FRACT; 

ERD; 

R  INTERFACE  REG  :  BEGIN 

I  MINUS  FORM  :=  ALPHA  FORM  *  (one  /  { (1  +  GAM'D  *  H  R  I) ) 

*  (two  7  GAMM  /HRI-l/RI); 

I  FORM  ;=  -one  +  ALPHA  FORM  *  Tone  /  H_R  IT  * 

((two  *  GAM4  /BRD-  (one-  GAMD-/  R  I) 

-  (HEAT  DEN  FRA£T"~ *  (FORM  /  two) )  ; 

I  PLUS  FORM  :=  XLRlffi  FORM  *  (SQR(GAMi)  / 

(H_E_I  *  (one  +  GAMD ) )  *  (two  /  H_R_I  +  one  /  R_I)  ; 

END; 

R_EXTERNAL_REG  :  BEGIN 

I_MINUS_PORM  :»  ALPHA  FORM  *  FRACT  -  (HEAT_DEN_FRACT 
*  temp  FORM  /  eight)  ; 

I_P0RM  :=  -one  *  ALPHA  FORM  *  (FRACT  +  (one  +  (two  /  H  R  I)) 

*  B  I)  -  (HEAT  DIN  FRACT  *  (three  *  temp  PORM/elght) ) ; 
I  PLUS_FOBM  :=  0.0; 

EHD;  ~ 

R  REG  END  :  BEGIN 

I  MINUS  FORM  :=  ALPHA  FORM  *  (one  /  H  R  I)  * 

(one  /  H“R  I  -  one  /  (two  T  R  I) ) ; 

I  FORM  :=  -one  *  ALPHK  Ft)RM  *  FRACT  -  (HEAT  DEN  FRACT  * 

W  b  nrml  [I  ^Ij]  /  two)  ; 

I_PLUS  FORM-  :=  ALPHA  FORM  *  one  /  H  R  I  *  (one  /  H  R  I  + 
one  r  (two  *  R_I) ) ; 

END; 

R  CENTR  END  :  BEGIN 

I_MINUS_PQRM  :=  0.0; 

I_FORM  :=  -one  *  ALPHA  FORM  *  FRACT  -  (HEAT  DEN  FRACT  * 

W  b  nrml  [I  VALT  /  tvro)  ; 

I_PLUS_FOrM  :=  ALPRA_PORM  *  FRACT; 

END; 

R  INTERFACE  END  :  BEGIN 

I  MINUS  FORM  :=  ALPHA  FORM  *  (one  /  ( (one  +  GAMD  *  H  R  I) )  * 

(two  *  GAM1  /  H  R  I  -  one  /  SI);- 
I  FORM  :=  -one  *  ALPHA  FORM  *  (one  7  H  R  I)  *  “ 

( (two  *  GANfl  /  H  R_I)  -  (one-“GSMD  /  R  I) 

-  (HEAT_DEN_FRAPT  *  W_b_nrml  [I_VAL]  7  two)  ; 

I  PLUS  FORM  :=  ALPHA  FORM  *  (SQR(GAMD  / 

(H_R_I  *  (one  +  GAMD ) )  *  (two  /  H_R_I  +  one  /  R_I) ; 

_  END; 

R  EXTERNAL  END  :  BEGIN 

I_M3NUS  FORM  :=  ALPHA  FORM  *  (two  /  SQR(H  R  I) )  - 

(HEAT_DEN  FRACT  *  W_bnrmT  [I  VAL]  /  eight)  ; 
IJPORM  :*=  -one  *  ALPHA  FORM  *  (FRACT  +  (one  +  two  /  H  R  I) 
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m 


*  B_i)  -  (  three  *  HEAT  DEN  FRACT  * 
W  b  nrml  [I  VAL]  /  eight)  ; 

I  PLUS  FORM-:*  0.0;  " 

E HD;  " 

END; 

END;  (*  rad_dir  *) 


PROCEDURE  S_VECTQR  <I_VAL,  J_VAL  :  INTEJGER;  VAR  S_VECT  :  GRID)  ; 

{ . calculates  s_vector  coefficient  in  preparation  for 

mesh  traversal  . } 


VAR 

A,  W  FORMULA,  T  FORMULA, 

HEATHDEN  FRACT,  ALPHA  FORM,  A  SQRD, 
HJR,~  N0CE_VAL  :  extended; 
teqp_u_v  : extended; 

BEX31N 

(*  estblish  node  type  *) 

IF  (J  VAL  =  M)  AND  (I  VAL  =  N)  THEN 
S  NODE  ;=  S  END  EXTERNAL 
ET.qff  ~  ~ 

IF  (J  VAL  =  M)  THEN 
S  NODE  :=  S_END 

else 

IF  (I  VAL  =  N)  THEN 
S  NODE  :=  S  EXTERNAL 
ktaV. 

S_NODE  :=  S  REG_CENTER  INT; 


(*  SJVECTOR  CALCULATIONS  -  PRELIMINARY  FORMULAS  *) 

A  SQRD  :  =  (BL  ALPHA  /  ALPHA  PROP  [I  VAL]  )  *  SQR(LEN/PHYS  RAD)  ; 
ALPHA  FORM  :=  ALPHA  PROP  [IlfflL]  /  BL  ALPHA; 

HEAT  DEN  FRACT  :=  (ORG  VALSTBBLOCD]  .DENSITY  * 

ORG  VALS  [BBLOCD]  .HEAT_CAP)  / 
(DENSrrY_AREAY[I_VAL]  *  HEAT_CAP  ARRAY  [I_VAL] )  ; 


H_r  :=  H_R__ARRAY  [I_VAL]  ; 

if  s_node  =  s_external  then 
begin 

tenpju  v  :=  0.0; 
w  formula  :=  0.0; 
end 
else 

VJPORMULA  :=  W_b_nrml  [I_VAL]  +  U_A_div; 

IF  J  VAL  =  M  THEN 

NODE  VAL  :=  T  A  NCDE[J  VAL  -1] 

RTCT.  — 

NODE_VAL  :=  T_A_NODE  [J_VAL] ; 

T_FORMULA  :=  NDDE_VAL  -  (TEMP_GRID  [I_VAL-1,  J_VAL]  /  eight)  ; 

(*  SJVECTOR  CALCULATIONS  -  FINAL  FORMULAS  *) 

CASE  S  NODE  OF 

S  REG  CENTER  INT  : 

S  VECT  [I  VAL,  J  VAL]  :  = 

HEAT  DEN  FRACT  *  (Q  ARRAY  nrml  [I  VAL]  + 

_  WJPORMULA  *  T_AjNCDE  [  JVAL]  +"U_V_dIv  *  TJVjNODE  [J  VAL] )  ; 

S  EXTERNAL  : 

S_VECT  [I  VAL, J  VAL]  := 

HEATDEN  FEaCT  *  (Q  ARRAY  nrml  [I  VAL]  + 

W  FORMULA  *  T~ FORMULA  +tenp  0  V  *  TT  V  NODE  [J  VAL]  - 
“TEMP  GRID  [T  VAL-1,  J  VALj  7  eight!)  ~+  ALPHA- FORM  * 

(one  +  two/  H  R)  *  B  I  T  0; 
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S  END 


S  VECT  [I  VAL, J  VAL]  := 

“  HEAT  DEN  FRACT  *  (Q  ARRAY  nrml  [I  VAL]  + 

vTb  nrml  [I  VAL]  *  T  A"NCDE  [J  val-1])  + 

<T  two  *  B_T_1  (I_VAL]~*  T_0)  7  (A_SQRD  *  H_Z) )  ; 

S  END  EXTERNAL  : 

~  "  S  VECT  [I  VAL,  J  VAL]  :  = 

"HEAT  DEfT  FRACT  *  (Q  ARRAY  nrml  [I  VAL]  + 

W  5  nrml  [I  VAL]  T  T_FOKMULA  )  + 

AEPHA_F0RR  *  (one  +  two/H  R)  *  B  I  *  T  0  + 

(two  *  B_i  1[I  VAL]  *  T  5)  /  (A"SQRD  7  H  Z,  ; 
END  (*  CASE  *) ;  ~ 

END;  (*  s  vector  *) 


PROCEDURE  ESTABLISH JX)EF_ARRAYS; 

{ . creates  arrays  of  axial  and  radial  coefficients  using 

arrays  of  records  . } 

VAR 

X,  J  :  INTEGER; 

BEGIN 

(*  establish  arrays  of  coefficients  in  the  radial  direction 
n.b.  --  each  I  value  references  a  record  containing 

three  coefficients  found  in  arrays  i  mat,  i  mat2, 
ORXG_A_z.  *)  ~ 

FOR  i  :=  n  downto  1  do 
FOR  j  :=  m  downto  2  do 

RAD  DIR  (I,  J,  ORIG_A  r  COEF  [I,  J]  .MINUS, 

~  ORIG  A  r  COEF  [T,  J]  .STD,  ORIG  A  r  COEF  [I,  J]  .PLUS)  ; 


FOR  I  :=  n  downTO  1  DO 
FOR  J  :=  m  downto  2  do 

AXIAL  DIR  (I,j,  ORIG  A  Z  COEF  [I, J]  .MINUS, 

0RTG_A_Z_00EF [I ,  JT.5TE,  ORIG_A_z_COEF [I,  J]  . PLUS)  ; 

END;  (*  establish_coef_arrays  *) 


PROCEDURE  TRAVERSE; 

{ . traverses  the  mesh  in  the  radial  and  axial  directions 

using  THCMAS '  ALGORITHM . } 

VAR 

INDEX,  X, start,  prev  total,  total  :  INTEGER; 

S  PTR,  S  TRV, 

CTPTR,  CTTRV, 

A- PTR,  A- TRV, 

B" PTR,  B"TRV, 

CT  PTR,  Cl  TRV, 

AI  PTR,  AITHV, 

BI  PTR,  BI~TRV, 

CXTPTR,  OTTRV, 

Ail" PTR,  AJ" TRV, 

BlT" PTR,  BJ~ TRV, 

D_PlR,  D_Tl?V,  NEWNQDE  :PTR_TYPE; 

S  V  COEF  :  GRID; 

ECMP,  OTMP3  :  CHAR; 


PROCEDURE  UPDATE_TISSUE_1  (VAR  TISSUE  :  GRID) ; 
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VAR 

I_OOUNT  :  INTEGER; 

BEGIN 

FOR  I  COUNT  :=  1  TO  N  DO 

Tissue  [i_count,i]  :=  30  /  normalizing  temp,- 


END; 

PROCEDURE  THCMAS_ALG(SYS_SIZE  :  INTEGER)  ; 


VAR 

I.ii.jj.Z  ••  INTEGER; 

P  PIR,  P  TRV, 

QT  PTR,  51  TRV, 

BK  Ql,  BK  5l  TRV, 

BlCP,  BKJPJIRV,  NEWNODE  :  PTRJTYPE  ; 

PROCEDURE  CAL  CPRIME  DPRIME; 

VAR  ~ 

I  :  INiEGER; 

FORM1  :  EXTENDED; 


BEGIN 

NEW (NEWNODE) ; 

NEWNODE"  .  INFO  :=  B  PTR".  INFO  /  A  PTR". INTO; 

NEWNODE". NEXT  :=  NIL; 

NEWNODE". BACK  ;=  NIL; 

P_PTR  :=  NEWNODE; 

P  TRV  :=  P  PTR; 

BK  P  :  =  NIL; 

NEW (NEWNODE) ; 

NEWNODE" . INTO  :=  D  PTR". INTO  /  A  PTR". INTO; 

NEWNODE". NEXT  :=  NIL; 

NEWNODE". BACK  :=  NIL; 

Ql  PTR  ;=  NEWNODE; 

QlTRV  :=  Ql  PTR; 

BK  Ql  :=  NIL“ 

AJTRV  :=  A  PTR". NEXT; 

C  TRV  :  =  Cf PTR". NEXT; 

B~TRV  :=  B_PTR" .NEXT; 

D  TRV  :=  D  PTR". NEXT; 

FOR  I  :=  2~T0  SYS  SIZE  DO 
BEGIN 

PORM1  :=  A  TRV". INTO  -  C  TKV" . INTO  *  P  TRV" . INTO; 

NEW  (NEWNODE) ; 

NEWNODE" . INTO  :=  B  TKV" . INTO  /  TORM1; 

NEWNODE". NEXT  :=  NIL; 

NEWNODE". BACK  :=  P  TRV; 

P  TRV". NEXT  :=  NEWJ30DE; 

PTRV  :=  P  TRV" .NEXT; 

NEW  (NEWNODE)  ; 

NEWNCDE". NEXT  :=  NIL; 

NEWNODE" . INTO  :=  (D  TRV".INro  -  C  IRV".3NTO  *  Ql  TRV" . INTO)  /  TORM1; 
NEWNODE" . BACK  :=  QlTKV; 

Ql  TRV"  .NEXT  ;=  NEWNODE; 

Ql- TRV  :=  Ql  TRV" .NEXT; 

A  TRV  :=  A  Tf?V"  .NEXT; 

CTTRV  :=  CT TRV". NEXT; 

B  TRV:=  B  TRV".NEXT; 

D_TRV:=  D~TRV" . NEXT ; 

END; 

BK  P  :=  P  TRV; 

BK_Q1  :=  5l_TKV; 

END; 

PROCEDURE  CALC  VALS; 

VAR 

I  :  INTEGER; 
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BEGIN 

NEW(NEWNODE) ; 

NEWNODE". INFO  :=  BK  Ql".INFO; 

NEWNODE". NEXT  :=  Nil; 

U  PTR  NEWNODE; 

BK  Q1  TRV  :=  BK  Ql".BACK; 

U  W :=  U  PTR;- 
BK_P_TRV  :=  BK_P" .BACK; 

(*  U  IS  BEING  FIT  .LED  VIA  FIRST  NODE  INSERTION  *) 

FOR  I  :=  SYS  SIZE-1  DOWNTO  1  DO 
BEGIN 

NEW  (NEWNODE)  • 

NEWNODEA . INFO  :=  BK  Q1  TRV".INPO  -  BK  P  TRV'MNPD  *  U  TRV" .  INPO; 
NEWNODE".NEXT  :=  U  PlRj 
U  PTO  :=  NEWNCDE;  “ 

BK  Q1  TRV  :=  BK  Q1  TRV" .BACK; 

BlO>_TRV  :=  BK_P_TEVa . BACK ; 
u_trv  :=  u_ptr; 

END; 

END; 

BEGIN  (*  THOMAS  ALG  *) 

MARK  (p); 

NEW  (U  PTR)  ; 

U  PTR'.INPO  :=  0.0; 
tT PTR A. NEXT  :=  NIL; 

NEW  (Q1  PTR) ; 

Q1  PTR'.INPO  :=  0.0; 

Q1  PTO".NEXT  :=  NIL; 

NEW  (P  PTR) ; 

P  PTR"7lNP0  :=  0.0; 

P  PTR A. NEXT  :=  NIL; 

T r TRV  :=  U  PTR; 

Ql  TRV  :=  Ql  PTR; 

P  5RV  ;=  P  PTR; 

FOR  I  :  =2  TO  SYS  SIZE  DO 
BEGIN 

NEW (NEWNODE) ; 

NEWNODEA . INFO  :=  0.0; 

NEWNODEA.NEXT  :=  NIL; 

U  TRV". NEXT  :=  NEWNODE; 

ITTRV  :=  U  TRV".NEXT; 

NEW  (NEWNODE)  ; 

NEWNODE" . INFO  ;=  0.0; 

NEWNODE". NEXT  :=  NIL; 

Ql  TRV". NEXT  :=  NEWNODE; 

Q1~TRV  :=  Ql  TRV". NEXT; 

NEW (NEWNODE) 7 
NEWNODE" . INFO  :=  0.0; 

NEWNODE" .NEXT  :=  NIL; 

P  TRV" .NEXT  :=  NEWNODE; 

PJTKV  :=  P  TRV". NEXT; 

END; 

CAL  CPRIME  DPRIME; 

CALCVALS; 

END; 


PROCEDURE  SEOCND_TRAVERSE_AXIAL ; 

{ . . .  traverses  the  mesh  in  the  axial  direction  moving  necessary 
coefficients  from  record  storage  structure  into  extended 
precision  arrays  for  access  by  Thomas'  Algorithm . } 


VAR 

I,J,j  find  :  INTEGER; 

I  MDrTEMP,  I  MAX  TEMP,  HALF  DELTA  :  extended; 
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INDEX  :  INTEGER; 

BEGIN 

MARK  (PI); 

S  PTR  :=  NIL; 

CTPTR  :=  NIL; 

A  PTO  :=  NIL; 

B~PTR  :*  NIL; 

Cl  PTR  ;=  NIL; 

AI"PTR  :=  NIL; 

BI— PTR  :=  NIL; 

D_PTR  :=  NIL; 

half  delta  :=  delta/2 ; 

TOTAE  :=  0; 

PREV  TOTAL  :=  1; 

FOR  T  :»  1  to  n  do 
BEGIN 
X  ;=  0; 

FOR  j  :=  2  to  m  do 
BEGIN 

X  ;=  X  +  1; 

(*  lead  coefficient  arrays  *) 

S_VECTOR  (I,  J,  S_V_COEF) ; 

NEW  (NEWNODE)  ; 

NEWNODE" .  INTO  :=  S  V  CJOEF  [I,J]  *  HALF  DELTA;  (*  S[X]*) 

NEWNODE" . NEXT  :=  NIL; 

IF  S  PTR  =  NIL  THEN 
BEGIN 

S  PTR  :=  NEWNODE; 

S^TRV  :=  SJPTR; 

END 

ELSE 

BEGIN 

S  TKV".NEXT  :=  NEWNODE; 

S_TRV  ;=  S_THV" . NEXT ; 

END; 

NEW  (NEWNODE) ; 

NEWNODE" . INFO  :=  -1  *  ORIG  A  Z  COEF  [I,  J]  .MINUS  *  HALF  DELTA; 
NEWNODE" .NEXT  :=  NIL; 

IF  C  PTR  =  NIL  THEN 
BEGIN 

C  PTR  :=  NEWNODE; 

CTTRV  :=  C_PTR; 

END 

ELSE 

BEGIN 

C  TRV" .NEXT  :=  NEWNODE; 

C_TRV  :=  C_TRV" . NEXT ; 

END; 

NEW  (NEWNODE)  * 

NEWNODE". INFO'  :=  -1  *  ORIG  A  z_COEF  [I,  J]  .STD  *  HALF  DELTA  + 
NEWNODE". NEXT  :=  NIL;  ~ 

IF  A  PTR  =  NIL  THEN 
BEGIN 

A  PTR  :=  NEWNODE; 

AJTRV  :=  A_PTR; 

END 

ft.se: 

BEGIN 

A  TRV". NEXT  :=  NEWNODE; 

A_TRV  :=  A_TRV" . NEXT ; 

END; 

NEW  (NEWNODE)  ; 

NEWNODE". INFO  :=  -1  *  ORIG_A_Z_OOEF  [I, J]. PLUS  *  HALFJDELTA ; 


(*C[X]*) 


1;  (*A[X]  *) 


(*  B[X]*) 
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NEWNCDE* .NEXT  :>  NIL; 

IF  B  PTR  »  NIL  THEN 
BEGIN 

B  PTO  :»  NEWNODE; 

STRV  B  PTR; 

END- 

ELSE 

BEGIN 

B  TRY* .NEXT  :=  NEWNODE; 

B— TRV  :=  B_TRV* . NEXT ; 

END; 

NEW  (NEWNCDE)  ; 

NEWNODE* . INFO  :=  ORIG  A  r  COEF  [I,  J]  .MINUS  *  HALF_DELTA;  (*  CI[X]*) 
NEWNODE*. NEXT  :=  NIL; 

IF  Cl  PTR  =  NIL  THEN 
BEGDT 

CI  PTR  :=  NEWNODE; 

CIjTRV  :=  CI_PTR; 

END 

tct.se: 

BEGIN 

Cl  TRV*  .NEXT  :=  NEWNODE; 

Cl  TRV  :=  Cl  TRV*. NEXT; 

END;- 

NEW  (NEWNODE)  ; 

NEWNODE* . INFO  :=  ORIG  A  r  COEF  [I,J]  .STD  *  HALFJDELTA  +  1;  (*  AI[X]*) 

NEWNODE*. NEXT  :=  NIL; 

IF  AI  PTR  =  NIL  THEN 
BEGIN- 

AI  PTR  :=  NEWNODE; 

AI'T’RV  :=  AI  PTR; 

END 

T7I.SK 

BEGIN 

AI  TRV*. NEXT  :=  NEWNODE; 

AI_TRV  ;=  AI_TRV* . NEXT ; 

END; 

NEW  (NEWNODE)  • 

NEWNODE* . INFO* :=  ORIG_A_r  COEF  tI,J] .PLUS  *  HALFJDELTA;  (*  BI [X]  *) 
NEWNODE*. NEXT  :=  NIL; 

IF  BI  PTR  =  NIL  THEN 
BBGIN- 

BI  PTR  ;=  NEWNODE; 

BIjTRV  ;=  BI_PTR; 

END 

ELSE 

BEGIN 

BI  TRV* .NEXT  :=  NEWNODE; 

BI— TRV  :=  BI  TRV*. NEXT; 

END; 

(*  calculate  constant  term  *) 

NEW  (NEWNODE); 

NEWNODE*. NEXT  :=  NIL; 

IF  D  PTR  =  NIL  THEN 
BEGIN 

D  PTR  :=  NEWNODE; 

DjTRV  :=  D_PTR; 

END 

ELSE 

BEGIN 

D  TRV*. NEXT  :=  NEWNODE; 

D— TRV  :=  D_TRV* .  NEXT  ; 

END; 

IF  I  =  1  THEN 

D  TRV* .  INFO  :=  AI  TRV*. INFO  *  HALF  ®ID  [I,J]  + 

BI  TRV* .^FO  *  half _grid  Tl+l,J]  +  S  TRV*. INFO 
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ETT.QTT 

IF  I  -  N  THEN 

D  TRV^.INFO  :*  Cl  TRV^.INPO  *  HALF  GRID  [1-1, J]  + 

~  AI  TRV.INPO- *  HALF  ®ID  [I,JT  +  S  TRV* .  INFO 
ELSE  ~  “ 

D  TRV^.INFO  Cl  THV\  INFO  *  half _grid  [i-1, j]  + 

~  AI  TKV* . INFO  ‘"HALF  OOD  [I, J]  +  bi  THV'.INPO 
*~half_grid  [i+l,jT  +  S_1R\T.INPO; 

(*  adjust  d[x]  for  j  ■  2  position  *) 

IF  J  -  2  THEN 

D_TRV*.INPO  D_TRVA  .  INFO  -  CJIHV*.  INFO  *  HALFJ3RID  [1,1]; 

(*  save  I,J  locations  to  facilitate  placing  of  new 
tenperatures  in  proper  locations  *) 

TOTAL  TOTAL  +  1; 

LOCATE  [TOTAL]. I  LOC  :=  I; 

LOCATE  [TOTAL]  .(TWC  :=  J; 

END;  (*  J  LOOP  *) 


THCMAS_ALG  (X)  ; 

(*  move  tenperatures  generated  by  THCMAS_ALG  to 
linear  array  locations  *) 

U  TRV  :=  U  PTR; 

FOR  INDEX  7=  PREV  TOTAL  TO  TOTAL  DO 
BEGIN 

LOCATE  [INDEX]  .TCM  TEMP  :=  U  TRV^.INFO; 

U_TRV  :=  U_TRV".NEXT; 

END; 

RELEASE  (P)  ; 

RELEASE  (PI); 

PREV  TOTAL  :=  PREV  TOTAL  +  X; 

END;  T*  I_LOOP  *)  ~ 

(*  move  temperatures  generated  by  IHCMAS_ALG  to 
proper  mesh  locations  *) 

FOR  INDEX  :=  1  TO  TOTAL  DO 

HALF_GRID  [LOCATE  [INDEX]  .I_LOC,  LOCATE  [INDEX]  .J_LOC]  :=  locate  [INDEX]  .tom_terp; 
DERIVE_CENTER_TEMP  (HALF_GRXD) ; 


END; 

PROCEDURE  FIRST_TRAVERSE_RADIAL; 

{ . . .traverses  the  mesh  in  the  radial  direction  moving  necessary 
coefficients  from  record  storage  structure  into  extended 
precision  arrays  for  access  by  Thomas '  Algorithm . } 


VAR 

I,J,J  FIND  :  INTEGER; 

J  MAX JTEMP, HALF JOELTA  :  extended; 
tT.tj  :  integer; 

INDEX  :  INTEGER; 

BEGIN 

S  PTR  :=  NIL; 

CTPTO  :=  NIL; 

A  PTR  :=  NIL; 

B  PTR  :=  NIL; 

CSJ  PTR  :=  NIL; 

Al~ PTR  :=  NIL; 

BJ~ PTR  :=  NIL; 

D_PTR  :=  NIL; 


-C36- 


MARK (PI) ; 

HALF  DELTA  :»  DELTA  /  2; 

PRE\TTCTAL  1; 

TOTAL  0; 

FOR  J  :«  2  TO  m  DO 
BEGIN 

X  0; 

FOR  I  1  to  n  DO 
BEGIN 

X  :«  X  +  1; 

(*  load  coefficient  arrays  *) 

S  VE CTOR  (I,  J,  S  V  COEF)  ; 
lyCW  (NEWNCDE)  * 

NEWNCDE'.  INFO '  :*  S  V  OOEF[I,J]  *  HALF  DELTA;  (*  S[X]  *) 
NEWNCDE'. NEXT  NTL7 

IF  S  PTR  ■  NIL  THEN 
BESIN 

S  PTR  :«  NEWNODE; 

S_TRV  :=  S_PTR; 

END 

PT.SF. 

BEGIN 

S  TRV'.NEXT  ;=  NEWNODE ; 

S~THV  :»  SJIRV'.NEXT; 

END; 


NEW  (NEWNCDE)  ; 

NEWNODE''. INFO  U  -1  *  ORIG  A  r  COEF  [I,J]  .MINUS  *  HALF  DELTA; 
NEWNODE'. NEXT  :=  NIL; 

IF  C  PTR  =  NIL  THEN 
BEGIN 

C  PTR  NEWNCDE; 

CTlHV  :=  C_PTR; 

END 

ELSE 

BEGIN 

C  TRV'.NEXT  :=  NEWNODE; 

C  TRV  :=  CJIRV'.NEXT; 

END; 

MEW  (MEWMODE)  ; 

NEWNCDE'. INFO'  :«=  -1  *  ORIG  A  r  COEF  [I,  J]  .STD  *  HALF  DELTA  +  1 
NEWNCDE'. NEXT  :=  NIL; 

IF  A  PTR  =  NIL  THEN 
BEGIN 

A  PTR  :=  NEWNCDE; 

A*~TKV  :=  A  PTR; 

END 

BEGIN 

A  TKV*  .NEXT  :=  NEWNODE; 

A  TRV  :=  A  TRV'.NEXT; 

END; 


ME!W  (NEWNCDE)  * 

NEWNCDE'. INFO'  :=  -1  *  ORIG  A  r  COEF  [I,  J]  .PLUS  *  HALF  DELTA; 
NEWNCDE'. NEXT  :=  NIL; 

IF  B  PTR  =  NIL  THEN 
BEGIN 

B  PTR  :=  NEWNCDE; 

BJJHV  :=  B_PTR; 

END 

ELSE 

BEGIN 

B  TRV'.NEXT  :=  NEWNODE; 

B~TRV  :=  B_TRV' . NEXT ; 

END; 


(*  C[X]*) 


;  (*  A[X)  *) 


(*  B  [X]  *) 
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NEW  (NEWNODE) ; 

NEWNODE" .  INTO  ORIG  A  Z  COEF  [I,J]  .MINUS  *  HALF  DELTA;  (*  C[J]*) 
NEWNODE" .NEXT  :«  NIL; 

IF  CJ  PTR  =  NIL  THEN 
BEGIN- 

CJ  PTR  :«  NEWNODE; 

CJTRV  ;=  CJ  PTR; 

END  ~ 

ELSE 

BEGIN 

CJ  TRV".NEXT  :=  NEWNODE; 

GTTRV  CJ_TKV"  . NEXT ; 

END; 


NEW  (NEWNCOE)  • 

NEWNODE" .  INFO'  ORIG  A  Z  COEF  II,  J]  .STD  *  HALF_DELTA  +  1;  (*  AJ[X]*) 

NEWNCOE ".NEXT  :=  NIL; 

IF  AJ  PTR  -  NIL  THEN 
BBGDT 

AJ  PTR  :=  NEWNODE; 

AJ“TRV  :=  AJ  PTR; 

END 

ELSE 

BEGIN 

AJ  TRV".NEXT  :=  NEWNODE; 

AJTRV  :=  AJ  TRV" . NEXT ; 

END; 


NEW  (NEWNCOE)  ; 

NEWNCDE" . INFO  :  =  ORIG_A  Z  COEF  [I, J], PLUS  *  HALF_DELTA ;  (*  BJ[X]*) 

NEWNODE". NEXT  ;=  NIL; 

IF  BJ  PTR  =  NIL  THEN 
BEGIN- 

BJ  PTR  :=  NEWNODE; 

BCTTRV  :=  BJ_FTR; 

END  ~ 

ELSE 

BEGIN 

BJ_TRV".NEXT  :=  NEWNODE; 

BJ_TRV  :=  BJ_TRV" .  NEXT  ; 

END; 


( ‘calculate  constant  term  *) 

NEW  (NEWNODE) ; 
newnode" . info  :=  0.0; 

NEWNODE". NEXT  :=  NIL; 

IF  D  PTR  =  NIL  THEN 
BEGIN 

D  PTR  :=  NEWNODE; 

D_TKV  :=  D_PTR; 

END 

ELSE 

BEGIN 

D  TRV" .NEXT  :=  NEWNODE ; 
DJIHV  :=  D_TRV" .  NEXT  ; 

END; 


IF  J  =  M  THEN 

D  TRV". INFO  :=  CJ  TRV".INFO  *  HALF  GRID[I,J-1]  + 
AJ  TRV". INFO  *  HALF_GRID  [I,J]  + 

S  iRTV"  .  INFO 


ELSE 

D_TRV" .  INFO  :=  CJ_TRV".INFO  *  HALF_GRID[I,  J-l]  + 
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AJ  TKV*.INFO  *  HALF  GRID  [I,J]  + 

BCTTRV*.INPO  *  halfjgrid  [i, j-t-1]  +  s_THV* . INFO; 


(*  save  I,J  locations  to  facilitate  placing  of  new 
tenperatures  in  proper  locations  *J 
TOTAL  :=  TOTAL  ^  1; 

LOCATE  [TOTAL] .1  LOC  I; 

LOCATE  [TOTAL]  .J-LOC  :=  J; 


END; 

THCMAS_AL3  (X)  ; 

(*  move  tenperatures  generated  by  7HCMAS_ALG  to 
linear  array  locations  *) 

U  TRV  ;=  U  PTR; 

FOR  INDEX  :=  PREV  TOTAL  TO  TOTAL  DO 
BEGIN 

LOCATE  [INDEX]  .TOM  TEMP  :=  U  TRV*.INPO; 

U_TRV  :=  U_TRV*  . NEXT ; 

END; 

RELEASE  (P)  ; 

RELEASE (PI) ; 

PREV  TOTAL  :=  PREV  TOTAL  +  X; 

END; 

(*  move  tenperatures  generated  by  THCMAS_ALG  to 
proper  mesh  locations  *) 

FDR  INDEX  :=  1  TO  TOTAL  DO 

HALFJ3RID  [LOCATE  [INDEX]  .  I _LOC, LOCATE  [INDEX]  .J_LOC]  :=  locate  [INDEX] .tom_tenp; 
DERIVE_CENTER_TSMP  (HALFjGRID) ; 


END; 


PROCEDURE  FILEJOUT 
VAR 


(STEP: REAL;  VAR  TISSUE : GRID; 

T_A_TEMP,  T_V_TEMP  :  J  ARRAY) ; 


INDEX  :  INTEGER ; 
ROW,  COL  :  integer; 
AVG  :  J  ARRAY; 
BEGIN 


(* 


(*  SENDS  OUTPUT  AT  IDENTIFIED  TIME  INTERVAL  TO  FILE  TEMP  OUT. DAT 
DATA  IS  PRINTED  AS  FORMATTED  REALS  IN  A  7:3  FIELD  USING 
A  COMA.  AS  THE  DELIMITER  *) 


IF  T  STEP  =  0.0  THEN 
BEGIN 

WRITELN  (CUT  DATA,  •  ORIGINAL  PROGRAM ')  ; 

WRITELN  (OOfDATA,  'ALL  VALUES  ARE  AT  THE  EXTERNAL  LAYER')  ; 

WRITE  (OUT  ETClA,  'STEP'  :10,  'J2'  :6,  'J6':12); 

WRITELN  (OUT  DATA, 'END1I' :  13,  'WMUS-6':10,  'WMUS-6 '  :  8,  'TA6  '  :5,  'TV6  '  :8)  ; 
END; 

WRITE  (OUT  DATA,  TSTEP:10:2,'  ,'); 

WRITE  (CXnTDATA,  HALF  GRID[N,  2]  *  NORMALIZING  TEMP:7 :3,  '  ,  ' )  ; 

WRITE  (OUT"DATA,  HALFJGRID [N, 6]  *NCRMALIZING  TEMP:7:3,  '  ,  ' )  ; 

WRITE  (CUTDATA,  HALFJGRID  [N,M]  *NCRMALIZING_TEMP:  7:3,  '  ,'); 

WRITE  (CUTDAIA,  W  ARKAY[6]  :7 :3,  '  ,',  W  ARR5Y[6]  :7:3,  '  ,'); 

WRITEIN  (COT  DATA,-!  A  NODE [6]  *  NORMALIZING  TEMP:7:3,  '  ,  ' 

,T  V  NODE  [6]  ‘NORMALIZING  TEMPT 7?3)  ; 
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*) 

(*  necessary  M  and  N  values  for  use  by  GRAPH  program  *) 
IF  T  STEP  =  0.0  THEN 
BEGIN 

WRITELN  (CUT  DATA,  M,  1  =M' )  ; 

WRITELN  (OOTTiATA.N,  '  =N' }  ; 

END; 


WRITELN  (COT  DATA,  STEP: 5:2)  ; 

WRITELN  (COT- DATA,  'T  A' )  ; 

TOR  INDEX  :=“1  TO  M-T  DO 
BEGIN 

WRITE  (COT  DATA,  T  A  TEMP  [INDEX]*  NORMALIZING  TEMP:7:3)  ; 

IF  INDEX  <>  M-l  THEN" 

WRITE  (COTJDAIA,  1  ,  ' )  ; 

END; 

WRITELN  (COT  DATA)  ; 

"  « 

WRITELN  (COT  DATA,  ’TV’); 

TOR  INDEX  :=  I  TO  M-I  DO 
BEGIN 

WRITE  (OOT  DATA,  T  V  TEMP  [INDEX]*  NORMALIZING  TEMP:7:3)  ; 

IF  INDEX  <>  M-l  THEN- 
WRITE  (OOTJDATA,  1  ,  * )  ; 

END; 

WRITELN  (COT_DATA)  ; 

WRITELN  (COT  DATA,  'TISSUE  TEMPERATURES’)  ; 

tor  row  :=  mxmro  1  DO 

BEGIN 

TOR  COL  :=  1  TO  M  DO 
BEGIN 

WRITE  (OOT  DATA,  TISSUE  [ROW,  COL]  *  NORMALIZING  TEMP  :  7:3); 
IF  COL  <>  R  THEN 

WRITE  (COT_DATA,  ’  ,  1 )  ; 

END; 

WRITELN  (0UT_DAIA); 

END; 

WRITEXN  (OOT  DATA,  ’EXTRAPOLATED  CENTER  LINE  TEMPERATURES')  ; 

FOR  COL  :=  1  TO  M  DO 

WRITE  (OOT  DATA,  TISSUE [0, COL] *  NORMALIZING  TEMP: 7: 3); 

WRITELN  (OOTJBATA)  ; 

WRITELN  (OOT  DATA)  ; 


END; 


BEGIN  (‘TRAVERSE  *) 

(*  initialize  step  counters  *) 

T  STEP  :=  0.0; 

AV  STEP  :=  0.0; 

HaIf_AV_STEP  :=  0.0; 

(*  pronpt  for  memory  dump  options  *) 


CLRSCR; 

WRITELN  (F,  '  INITIAL  VALUES  ' )  ; 

WRITEIN  (F)  ; 

DUMP  HALF  ARRAYS; 

DUMP2r_AREAY  (T_STEP,  TEMP_GRID)  ; 

(*  INITIALIZE  GRID  *) 

HALF  GRID  :=  TEMP  GRID; 

HALEr"AV  STEP  :=  HALF  AV  STEP  +  0.5; 

FILETCOT  (T  STEP,  HALF  GRID,  T  A  NODE,  T  V  NODE)  ; 
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CALCJIAJTVHALF  (T_STEP) ; 


REPEAT 

RET  EASE  (HEAPORG)  ; 

T_STEP  :=  T_STEP  +  0.5; 

CALC_TA_TV_HALF  (T_STEP) ; 

FIRST_TRAVERSE__RADIAL  ; 

RET  EASE  (HEAPORG); 

IF  T  STEP  =0.5  THEN 
BEGIN 

DUMP  HALF  ARRAYS; 

DUMP-T_ARRAY  (T_STEP,  HALFjGRID) ; 
END; 

T_STEP  :=  T_STEP  +  0.5; 

SECOND  TRAVERSE  AXIAL; 

RET  EASE  (HEAPORG)  ; 


TEMP  GRID  :=  HALF  GRID; 

HALF- AV_STEP  :=  HSLF_AV_STEP  +  1; 
AV_STEP  :=  AV_STEP  +  1; 

T_a  :=  T_a  bar; 

T  v  :=  T_v~bar; 

oTd  t_v_no3e  :=  t_v_node; 

old~t  a  node  :=  t  a  node; 


(*  print  time  step  counter  *) 

CLRSCR; 

GOTOXY  (40,12); 

WRITELN  (trunc  (t_step) )  ; 

IF  (T  STEP/INTERVAL  =  TRUNC (T  STEP  /  interval)  *  1.0) 
BESlN 
CLRSCR; 

DUMP_HALF  ARRAYS; 

DUMP_T  ARRAY  (T  STEP,  HALF  GRID)  ; 

FILEJXJT  (T_STEP,  HALF_GRID,T_A_NQDE,  T_V_NODE)  ; 
end; 


(*  DUMP  DATA  FOR  GRAPHICS  TEST  PROGRAM  *) 

(*  IF  T_STEP  =  1000  THEN 

filejout  (t_step,  half_grid,  t_a_node,  t_v_node) ;  *) 

(♦RETURN  TISSUE  TEMP  AT  J=1  TO  30 
IF  T  STEP  =  100  THEN 

Update  tissue  i  (half  grid) 

*)  “ 


UNTIL  LCV_R  =  T_STEP; 
END;  (*  TRAVERSE  *) 


BEGIN  (*  TISSUE  *) 

ESTABLISH_COEF_ARRAYS ; 
TRAVERSE; 

END;  ( ‘TISSUE*) 

end. 


THEN 
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4.  PHYS_SP.PAS 

{  ANATOMICAL  INPOT  MODULE  } 


UNIT  PHYS_sp; 


INTERFACE 
USES  var  sp; 

PROCEDURE  ANATOMICAL  INPOT  (VAR  DIA,  LEN  :  extended)  ; 

PROCEDURE  NUM_DIVI S ICNS  (VAR  CORE  SUBDIV,  MUSCLE  SUBDIV, 

EAT  SUBDIV,  SKIN  SUBDIV, ~N,  M  :  INTEGER)  ; 

PROCEDURE  WIDIHJCALCUIATTCN5  (VAR  H  CORE,  H  MUS,  H  EAT,  H  SKIN, 

C5RE_MUS  SAMMA,  HUS  EAT  GAMMA, 

EAT  SKIN  GAMMA  :  extended; 

DIA,  LEU  :  REAL; 

CORE  DIV,  MUS  DIV,  EAT  DIV, 

SKHTDIV  :  INTEGER)  ;  ” 

PROCEDURE  EST_INTEREACES  (VAR  BOUND  SET  ;  BOUNDS; 

VAR  CORE  MUS_BNDRY,  MUS  FAT  BNDRY, 

FAT  SKUsT BNDRY,  SKIN  SURFACE:  INTEGER; 

CORE  SUBDIV,  "MUS  SUBDIV,  ~FAT  SUBDIV, 

Skin  subdiv  :  integer)  ,•  ~ 

PROCEDURE  CREATE_ARRAYS  (VAR  H  ARRAY,  DISTANCE_ARRAY  :  I_ARRAY; 

CORE  MUS  BOUND,  MUS  FAT  BOUND, FAT  SKIN  BOUND: INTEGER; 
H  r  CORE,  H  r  MUS,~H  r  EAT,  H  r  SKIN  :  REAL)  ; 

IMPLEMENTATION 

USES 

CRT,  DOS,  PRINTER  ; 

(*=============================  ANATOMICAL  PROCEDURES  ======================*) 


PROCEDURE  ANATOMICAL__INPOT  (VAR  DIA,  LEN  :  extended)  ; 

{ . provides  for  keyboard  input  of  length  and  diameter  in  can. 

exports  length  and  diameter  in  meters  . } 

VAR 

RESPONSE  :  CHAR; 

BEGIN 

dia  :=  GC  DIA; 
len  :=  GC  LEN; 

{ 

CURSCR;  _  _ 

WRITE  (’ENTER  THE  DIAMETER  (cm)  :  ’)  ; 

READIN  (DIA)  ; 

WRITELN; 

WRITE  (’ENTER  THE  LENGTH  (cm)  :  ’); 

READIN  (LEN)  ; 

REPEAT 

CLRSCR; 

WRITELN  (’EXIT  /  VERIFY  SCREEN’); 

WPTTFTJJ- 

WRITELN*  ( ’  1)  DIAMETER  =  ’  :20,  DIA  :  15:7  ,  ’  cm’); 

WRITELN  (’2)  LENGTH  =  ’  :20,  LEN  :  15:7  ,  ’  cm’); 

WRITELN; 

WRITELN  (’ENTER  THE  NUMBER  OF  THE  ITEM  TO  BE  CHANGED’); 

WRITE  {’--- . ENTER  V  TO  VERIFY  AND  CONTINUE . >  ’)  ; 

READIN  (RESPONSE)  ; 

IF  (RESPONSE  =  ’1’)  OR  (RESPONSE  =  ’2’)  THEN 
BEGIN 
WRITELN; 

CASE  RESPONSE  OF 

'1'  :  BEGIN  _ 

WRITE  ('ENTER  THE  NEW  DIAMETER  (cm)  :  ')  ; 

READIN  (DIA)  ; 

WRITELN; 

END; 

'2'  :  BEGIN 
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WRITE  ( 'ENTER  THE  NEW  LENGIH  (cm)  :  '); 
READIN  (LEN) ; 

WRITELN; 

END; 

END;  (*  CASE  *) 

END; 

UNTIL  ((RESPONSE  =  'V')  OR  (RESPONSE  =  'v')); 


(*  CONVERSION  TO  METIERS  *) 

LEN  :=  LEN  /  100; 

DIAM  :=  DLAM  /  100; 

RAD  :=  DIA/2; 

PHYS  RAD  :=  RAD; 


END;  (*  ANATCMICAL_INPOT  *) 


PROCEDURE  NUM  DIVISIONS  (VAR  03RE  SUBDIV,  MUSCLE  SUBDIV, 

FATJSUBDlV,  SKXN_SUBDlV ,  N,  M  :  INTEGER)  ; 


{ . provides  for  user  entry  for  the  number 

of  divisions  in  each  organ . 


} 


VAR 

RESPONSE  ;  CHAR; 

BEGIN 

axial_divs  :=  GC_AXIAL_DIVS; 
m  :=  axial_divs  +  1; 
core  subdiv  :  =  GC  CORE  DIVS; 
muscXe  subdiv  :=  SC  MI®  DIVS; 
fat  subdiv  :  =  GC  EAT  DIVS; 
skin_subdiv  :=  GCjSKlN  DIVS; 

n  :=  core  subdiv  +  muscle_subdiv  +  fat_subdiv  +  skinjsubdiv  +  1; 

N  :=  9999; 

REPEAT 

CLRSCR; 

WP I  ’lTi'l  • 

WRITELN*  ( 'MAXIMUM  AXIAL  DIVISIONS  =  ' ,  J_MAX-1)  ; 

• 

WRITE  ('ENTER  THE  NUMBER  OF  AXIAL  DIVISIONS  DESIRED  :  ')  ; 
READIN  (AXIAL  DIVS)  ; 

WPTTRfJtt- 

IF  (AXIAL  DIVS  >  J  MAX-1)  OR  (AXIAL  DIVS  <=0)  THEN 
BEGIN 

WRITELN  ('AXIAL  DIVISIONS  EXCEED  PROGRAM  PARAMETERS')  ; 
WRITELN  ('AXIAL  INPUT  SEQUENCE  WILL  BEGIN  AGAIN’)  ; 
PAUSEjCRT; 

END; 

IF  (AXIAL  DIVS  >=1)  AND  (AXIAL  DIVS  <=  J  MAX-1)  THEN 
BEGIN  ~ 

CLRSCR; 

WRITELN; 

WRITE IN  ('VERIFY  /  EDIT  SCREEN'); 

WPT^'lfii  JJ  • 

WRITELN'  ( '  AXIAL  DIVISIONS  =  ' ,  AXIAL  DIVS)  ; 

M  :=  AXIAL_DIVS  +  1; 

WRITELN; 

WRITEIN  ('  M  =  1  ,M) ; 

WPTTVTJJ* 

WRITELN'  ( 'ENTER  V  TO  VERIFY  AND  CONTINUE  OR ')  ; 

WRITE  ( '  PRESS  E  TO  CHANGE  AXIAL  DIVISIONS  DATA  ' )  ; 

READIN  (RESPONSE); 

WRITEIN; 

CLRSCR; 

END; 

UNTIL  ( (RESPONSE  *  'V' )  OR  (RESPONSE  =  'v' ) )  ; 

RESPONSE  :=  '  '; 
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WHILE  N  >  I  MAX  DO 
BEGIN 

WRITELN  { 'ENTER  THE  NUMBER  OF  SUBDIVISIONS  FOR  THE  FOLLOWING  LAYERS  - ' )  ; 
WRITEIN  ('<  MAXIMUM  TOTAL  SUBDIVISIONS  =  ':50,  I_MAX-1,  '  )'); 

WRITEIN; 

WRITE  {'CORE  :  25); 

READIN  (CORE_SUBDIV) ; 

WRITELN; 

WRITE  ('MUSCLE  :  ’  :  25); 

READIN  (MUSCLE_SUBDIV)  ; 

WRITELN; 

WRITE  ( 1  EAT  :  *  ;25) ; 

READIN  (FAT  SUBDIV)  ; 

WRITELN; 

WRITE  ('SKIN  :  '  :25) ; 

READIN  (SKIN_SUBDIV)  ; 

WRITELN; 

REPEAT 

CLRSCR; 

WRITELN; 

WRITELN  ('VERIFY  /  EDIT  SCREEN'); 

WRTTHf  N  • 

WRITEIN'  ( 'SUBDIVISION  VALUES  ENTERED  ' )  ; 

WRITELN; 

WRITEIN  Cl)  CORE  :  '  :40,  CORE  SUBDIV); 

WRITELN  C2)  MUSCLE  :  '  :40,  MUSdE  SUBDIV); 

WRITEIN  C3)  FAT  :  '  :40,  FAT  SUBDIV)  ; 

WRITELN  ('4)  SKIN  :  ’  :40,  SKIN  SUBDIV)  ; 

WRITEIN; 

WRITEIN  ( 'ENTER  THE  NUMBER  OF  THE  LAYER  TO  BE  CHANGED' )  ; 

WRITE  (' . ENTER  V  TO  VERIFY  AND  CONTINUE . . >  ')  ; 

READIN  (RESPONSE)  ; 

VJRT'lv.i  N • 

IF  (RESPONSE  >=  '1')  AND  (RESPONSE  <='4')  THEN 
BEGIN 

WRITEIN; 

WRITE  ( 'ENTER  THE  NEW  VALUE  FOR  THE  ' )  ; 

CASE  RESPONSE  OF 
'1'  :  BEGIN 

WRITE  ('CORE  :  '); 

READIN  (CORE  SUBDIV) ; 

END; 

'2'  :  BEGIN 

WRITE  ('MUSCLE  :  ') ; 

READIN  (MUSCLE_SUBDIV)  ; 

END; 

'3'  :  BEGIN 

WRITE  ( ' FAT  ;  ' )  ; 

READIN  (FAT  SUBDIV)  ; 

END; 

'4'  ;  BEGIN 

WRITE  ( '  SKIN  :  ' )  ; 

READIN  (SKIN_SUBDIV)  ; 

END; 

END;  {*  CASE  *) 

WRITEIN; 

END;  (*  IF  *) 

UNTIL  ((RESPONSE  =  'V')  OR  (RESPONSE  =  'v')); 

N  ;=  CORE  SUBDIV  +  MUSCLE  SUBDIV  +  FAT  SUBDIV  +  SKIN  SUBDIV  +  1; 

IF  N  >  I  FJAX-1  THEN 
BEGIN  “ 

WRITEIN  ('DIVISION  TOTAL  EXCEEDS  PROGRAM  PARAMETERS’)  ; 

WRITEIN  ( 'MAXIMUM  TOTAL  DIVISIONS  CANNOT  EXCEED  I_MAX-1)  ; 

WRITELN; 

WRITEIN  ('DATA  ENTRY  SEQUENCE  WILL  START  AGAIN')  ; 

PAUSE_CRT; 

END; 

END;  (*  WHILE  *) 

(*  ESTABLISH  H_z  AND  NORMALIZE  *) 

PHYS_H_Z  :=  LEN  /  AXIAL_DIVS; 
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H_Z  : «PHYS_H_Z  /  LEN; 

END;  (*  num_di vi s i ons  *) 

PROCEDURE  WIDTH  CALCUIATICNS  (VAR  H  CORE,  H  MUS,  H  EAT,  H  SKIN,  CORE  MUS  GAM4A, 

MUS  EAT  G3&MA,  -  -  - 

EAT- SKIN  GAMMA  :  extended; 

DIA,  LESJ  :  REAL;  CORE  DIV,  MUS  DIV,  EAT  DIV, 
SKXN_DIV  :  INTEGER) ;  “ 

{ . inports  the  location  of  organ  boundaries,  length,  diameter 

calculates  width  of  each  organ 

exports  organ  width,  H-  values,  and  garrna  values . } 

BEGIN  (*  width_calculations  *) 


(*  establish  sanple  boundary  *) 
core  mus  bound  :=  core_ratio; 
musJEat_5ound  :=  muscle_ratio; 
fat_skin  bound  :=  fat_ratio  ; 
skin_surface_bound  :  =  skin_ratio  ,- 

(*  DETERMINE  LAYER  WIDTHS  *) 

CORE  WIDTH  :=  CORE  MUS  BOUND  ; 

MUS_WlDTH  :=  (MUSjRAT  BOUND  -  CORE  MUS  BOUND  )  ; 

EAT  WIDTH  :=  (EAT  SKIN  BOUND  -  MUS  EAT~BOUND)  ; 

SKIN  WIDTH  :=  (SKlN  SURFACE  BOUND  ~  EAT  SKIN  BOUND)  ; 


(*  DETERMINE  H-  VALUES  *) 

H  CORE  :=  (CORE  WIDTH  /  (CORE  DIV  +  0.5)); 

iTMus  :=  (mus  Width  /  mus  diV)  ; 

H- FAT  :=  (EAT  WIDTH  /  EAT- DIV)  ; 

ITSKIN  :=  (SKIN_WIDTH  /  SKlN_DIV)  ; 

(*  DETERMINE  GAW4A  VALUES  *) 

CORE  MUS  GAIWA  :=  H  CORE  /  H  MUS; 

mus  Rat  Gamma  :=  h  Mus  /  h  bet  ; 

EAT_SKlNjGAMMA  :=  H_EAT  /  H_SKIN; 

END;  (*  width_calculatians  *) 

PROCEDURE  EST  INTERFACES  (VAR  BOUND  SET  :  BOUNDS; 

VAR  CORE  MUS  BNDRY,  MUS  EAT  BNDRY, 

EAT  BKQrBNDRY,  SKIN  SURFACE:  INTEGER; 

CORE  SUBDIV,-MUS  SUBDIV,-EAT  SUBDIV, 

BKIN_SUBDIV  :  INTEGER)  ; 

{ . establishes  a  SET  of  I  values  representing  the  location  of 

organ  interfaces . } 

BEGIN 

CORE  MUS  BNDRY  :=  CORE  SUBDIV  +  1; 

MUS  TAT  BNDRY  :=  MUS  SUBDIV  +  CORE  MUS  BNDRY; 

EAT- SKIN  BNDRY  :=  MUS  EAT  BNDRY  +  RAT  SUBDIV; 

SKIN  SURRACE  :=  FAT  SKIN  BNDRY  +  SKUTSUBDIV; 

B0UNB_SET  :=  [CaRE_MUS_BNDRY ,  MUS_EAT^BNDRY,  EAT_SKIN_BNDRY]  ; 

END; 


PROCEDURE  CEEATE_ARRAYS  (VAR  H  ARRAY, DISTANCE  ARRAY  :  I  ARRAY; 

CORE  MUS_BCUND,MUS  EAT  BOUND,  PET  SKIN  BOUND :  INTEGER ; 
H_rJSDRE,  H_r_MUS,~H_r“FAT,  H_r_SKIN  T  REAL)  ; 

{ . creates  a  set  of  I  values  for  each  organ  and  establishes 

a  series  of  arrays  containing  R  i  and  H_r  values 

indexed  to  the  I  locations . } 
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VAR 

INDEX  :  INTEGER; 

H_TOTAL  :  REAL; 

BEGIN  (*  create  arrays  *) 

OGRE  SET  :»  [1. . (CORE  MUS  BOUND)] ; 

MUS  SET  :=  [CORE  MUS  BOUND  +  1.  .MUS  FAT  BOUND]  ; 

FAT  SET  :=  [MUS  PAT  BOUND  +  1. .FAT  SKIN  BOUND] ; 

SKIN  SET  [FAT  SKIN  BOUND  +  X.  .NT; 

FOR  INDEX  :=  1  TO  N  DC 
BEGIN 

IF  INDEX  IN  CORE  SET  THEN 
H  ARRAY  [INDEX!  :=  H  r  CORE; 

IF  INDEX  IN  MUS  SET  THEN 
H  ARRAY  [INDEX]  ;=  H  r  MUS; 

IF  INDEX  IN  FAT  SET  THEN 
H  ARRAY [INDEX!  :=  H  r  FAT; 

IF  INDEX  IN  SKIN  SET  THEN 
H_ARRAY  [INDEX!  :=  H_r_SKIN; 

END; 

H  ARRAY [0]  :=  0.0; 

DISTANCE  ARRAY  [0]  :=  0.0; 

H_ARRAY[T]  :=  H_ARRAY [2] /2 ; 

(*  establish  R_i  and  physical_r_i  arrays  *) 

FOR  INDEX  :=  1  TO  N  DO 

FHYS_H_r_ARRAY  [INDEX]  ;=  H_ARRAY  [INDEX]  *  PHYS_RAD; 

DISTANCE  ARRAY  [l]  ;=  H  ARRAY [1]  ; 

PHYS  R  i  [1]  :=  EHYS_H_r_ARRAY  [1] ; 
phys  r  i  To!  ;=  0; 

~Ft3R  INDEX  :=  2  TO  N  DO 
BEGIN 

DISTANCE_ARRAY  [INDEX]  :=  DISTANCE  ARRAY  [ INDEX- 1]  +  H  ARRAY  [INDEX]  ; 
FHYS_R_i  [INDEX]  :=  DISTANCE_ARRAY [INDEX]  *  PHYS_RAD7 
END; 

H  TOTAL  :=  0.0; 

FOR  INDEX  :=  1  TO  N  DO 

HJIOTAL  :=  HJTOTAL  +  H_ARRAY [INDEX]  ; 

END;  (*  createarrays  *) 

END. 
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APPENDIX  D  -  DISTRIBUTION  LIST 

2  Copies  to: 

Defense  Technical  Information  Center 
ATTN:  DTIC-DDA 
Alexandria,  VA  22304-6145 

Office  of  the  Assistant  Secretary  of  Defense  (Hlth  Affairs) 
ATTN:  Medical  Readiness 
Washington,  DC  20301-1200 

Commander 

U.S.  Army  Medical  Research  and  Development  Command 
ATTN:  SGRD-PLC 
Fort  Detrick 

Frederick,  MD  21702-5012 
Commander 

U.S.  Army  Medical  Research  and  Development  Command 
ATTN:  SGRD-PLE 
Fort  Detrick 

Frederick,  MD  21702-5012 
Commandant 

Army  Medical  Department  Center  and  School 
ATTN:  HSMC-FR,  Bldg.  2840 
Fort  Sam  Houston,  TX  78236 

1  Copy  to: 

Joint  Chiefs  of  Staff 

Medical  Plans  and  Operations  Division 

Deputy  Director  for  Medical  Readiness 

ATTN:  RAD  Smyth 

Pentagon,  Washington,  DC  20310 

HQDA 

Office  of  the  Surgeon  General 
Preventive  Medicine  Consultant 
ATTN:  SGPS-PSP 
5109  Leesburg  Pike 
Falls  Church,  VA  22041-3258 

HQDA 
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Assistant  Secretary  of  the  Army  for  Research,  Development  and  Acquisition 

ATTN:  SARD-TM 

Pentagon,  Washington,  DC  20310 

HQDA 

Office  of  the  Surgeon  General 
ATTN:  DASG-ZA 
5109  Leesburg  Pike 
Falls  Church,  VA  220^1-3258 

HQDA 

Office  of  the  Surgeon  General 
ATTN:  DASG-DB 
5109  Leesburg  Pike 
Falls  Church,  VA  22041-3258 

HQDA 

Office  of  the  Surgeon  General 
Assistant  Surgeon  General 
ATTN:  DASG-RDZ/Executive  Assistant 
Room  3E368,  The  Pentagon 
Washington,  DC  20310-2300 

HQDA 

Office  of  the  Surgeon  General 
ATTN:  DASG-MS 
5109  Leesburg  Pike 
Falls  Church,  VA  22041-3258 

Uniformed  Services  University  of  the  Health  Sciences 
Dean,  School  of  Medicine 
4301  Jones  Bridge  Road 
Bethesda,  MD  20814-4799 

Uniformed  Services  University  of  the  Health  Sciences 
ATTN:  Department  of  Military  and  Emergency  Medicine 
4301  Jones  Bridge  Road 
Bethesda,  MD  20814-4799 
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Commandant 

Army  Medical  Department  Center  &  School 
ATTN:  Chief  Librarian  Stimson  Library 
Bldg  2840,  Room  1 06 
Fort  Sam  Houston,  TX  78234-6100 

Commandant 

Army  Medical  Department  Center  &  School 
ATTN:  Director  of  Combat  Development 
Fort  Sam  Houston,  TX  78234-6100 

Commander 

U.S.  Army  Aeromedical  Research  Laboratory 

ATTN:  SGRD-UAX-SI 

Fort  Rucker,  AL  36362-5292 

Commander 

U.S.  Army  Medical  Research  Institute  of  Chemical  Defense 
ATTN:  SGRD-UVZ 

Aberdeen  Proving  Ground,  MD  21010-5425 
Commander 

U.S.  Army  Medical  Materiel  Development  Activity 
ATTN:  SGRD-UMZ 
Fort  Detrick 

Frederick,  MD  21702-5009 
Commander 

U.S.  Army  Institute  of  Surgical  Research 

ATTN:  SGRD-USZ 

Fort  Sam  Houston,  TX  78234-5012 

Commander 

U.S.  Army  Medical  Research  Institute  of  Infectious  Diseases 
ATTN:  SGRD-UIZ-A 
Fort  Detrick 

Frederick,  MD  21702-5011 
Director 

Walter  Reed  Army  Institute  of  Research 

ATTN:  SGRD-UWZ-C  (Director  for  Research  Management) 

Washington,  DC  20307-5100 
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Commander 

U.S.  Army  Natick  Research,  Development  &  Engineering  Center 
ATTN:  SATNC-Z 
Natick,  MA  01760-5000 

Commander 

U.S.  Army  Natick  Research,  Development  &  Engineering  Center 
ATTN:  SATNC-T 
Natick,  MA  01760-5002 

Commander 

U.S.  Army  Natick  Research,  Development  &  Engineering  Center 
ATTN:  SATNC-MIL 
Natick,  MA  01760-5040 

Commander 

U.S.  Army  Research  Institute  for  Behavioral  Sciences 
5001  Eisenhower  Avenue 
Alexandria,  VA  22333-5600 

Commander 

U.S.  Army  Training  and  Doctrine  Command 
Office  of  the  Surgeon 
ATTN:  ATMD 

Fort  Monroe,  VA  23651-5000 
Commander 

U.S.  Army  Environmental  Hygiene  Agency 
Aberdeen  Proving  Ground,  MD  21010-5422 

Director,  Biological  Sciences  Division 
Office  of  Naval  Research  -  Code  141 
800  N.  Quincy  Street 
Arlington,  VA  22217 

Commanding  Officer 

Naval  Medical  Research  &  Development  Command 
NNMC/Bldg  1 

Bethesda,  MD  20889-5044 
Commanding  Officer 

U.S.  Navy  Clothing  &  Textile  Research  Facility 

P.O.  Box  59 

Natick,  MA  01760-0001 


-  D4- 


Commanding  Officer 

Navy  Environmental  Health  Center 

2510  Walmer  Avenue 
Norfolk,  VA  23513-2617 

Commanding  Officer 

Naval  Aerospace  Medical  Institute  (Code  32) 
Naval  Air  Station 
Pensacola,  FL  32508-5600 

Commanding  Officer 

Naval  Medical  Research  Institute 

Bethesda,  MD  20889 

Commanding  Officer 
Naval  Health  Research  Center 
P.O.  Box  85122 
San  Diego,  CA  92138-9174 

Commander 

Armstrong  Medical  Research  Laboratory 
Wright-Patterson  Air  Force  Base,  OH  45433 

Strughold  Aeromedical  Library 
Document  Services  Section 

251 1  Kennedy  Circle 
Brooks  AFB,  TX  78235-5122 

Commander 

US  Air  Force  School  of  Aerospace  Medicine 
Brooks  Air  Force  Base,  TX  78235-5000 

Director 

Human  Research  &  Engineering 
US  Army  Research  Laboratory 
Aberdeen  Proving  Ground,  MD  21005-5001 
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